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Abstract 

A  higher-order  shell  theory  is  used  to  analyze  compressive  and  tensile  loads 
on  a  graphite/epoxy  laminated  cylinder  containing  a  square  cutout.  The  Hashin 
failure  criterion  is  used  to  determine  failure  in  the  fiber,  matrix,  or  lamination. 
Once  failure  occurs,  the  appropriate  stiffness  terms  are  reduced.  This  failure  causes 
a  redistribution  of  stress,  leading  to  further  failure.  In  order  to  account  for  the  loss  of 
residual  strength  due  to  cyclic  loading,  the  stiffness  matrix  is  further  reduced  at  each 
new  increment  of  load  or  displacement.  The  objective  is  not  to  determine  the  S-N 
curve  for  the  material,  but  rather  to  determine  the  damage,  displacement,  and  stress 
distribution  in  a  complex  configuration  under  fatigue  loading  using  a  progressive 
failure  and  stiffness  reduction  approach.  The  failure  progression  for  static  loads  is 
determined  and  used  as  an  indicator  for  the  cyclic  loads.  It  is  shown  that  as  the 
stiffness  decreases,  the  global  displacements  increase,  resulting  in  more  failure  until 
the  entire  panel  fails. 


FINITE  ELEMENT  ANALYSIS 
OF  A  COMPOSITE 
CYLINDRICAL  SHELL  WITH  A 
CUTOUT  UNDER  FATIGUE 

LOADING 


I.  Introduction 

1.1  Motivation 

A  common  problem  among  pressurized  cylinders,  such  as  aircraft  fuselages,  are 
fatigue  cracks  initiating  from  cutouts.  These  cutouts  are  often  rectangular  in  shape, 
such  as  doors  and  windows  in  a  fuselage,  resulting  in  large  stress  concentrations  at 
the  corners.  A  better  understanding  of  the  stress  and  damage  due  to  the  presence 
of  cutouts  will  hopefully  lead  to  improvements  in  composite  structural  design,  and 
this  research  specifically  targets  cylindrical  shells  due  to  their  common  structural 
applications. 

Fiber-matrix  composite  materials  are  being  increasingly  used  in  structural  ap¬ 
plications,  including  ones  involving  the  above-mentioned  geometry.  Their  tailorable 
properties  result  in  a  high  strength-to-weight,  ratio,  which  is  a  performance-driver 
in  many  applications.  They  do  not  often  behave  as  homogeneous  materials  do,  re¬ 
quiring  new  techniques  and  research.  For  reference,  the  term  “composite  in  this 
thesis  is  meant  to  describe  fiber-matrix  laminates,  where  parallel  continuous  fibers 
are  surrounded  by  a  matrix.  This  defines  a  unidirectional  laminae,  and  these  may  be 
stacked  in  different  fiber  orientations  to  develop  a  multidirectional  laminate.  Unless 
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otherwise  mentioned,  the  matrix  is  assumed  to  be  constructed  of  resin  or  epoxy  and 
not  metal. 

One  of  the  first  steps  in  design  is  to  account  for  the  static  loading.  Many  com¬ 
ponents,  however,  fail  in  fatigue  rather  than  static  loading.  Techniques  to  analyze 
fatigue  for  isotropic  materials  are  not  always  applicable  to  composite  materials,  and 
new  methods  must  be  explored. 

1.2  Background 

1.2.1  Shells.  Classical  plate  theory  falls  short  in  capturing  the  state 
of  stress  and  strain  in  a  curved  shell,  requiring  the  development  of  shell  theory. 
Shell  theory  attempts  to  simplify  this  analysis  by  modeling  thin  structures  as  two- 
dimensional  problems.  This  is  done  by  assuming  the  in-plane  loads  are  dominant, 
which  allows  the  use  of  a  datum  surface  to  characterize  out-of-plane  loads.  In  this 
way,  the  three-dimensional  problem  is  represented  two-dimensionally.  Because  of  the 
complicated  nature  of  curved  shells,  several  different  theories  have  been  proposed. 

One  of  the  first  shell  theories  was  the  Love  (1)  theory,  which  applied  Kirchhoff 
type  assumptions  to  a  shell.  This  is  a  classical  linear  approach.  Donnell  (2)  applied 
the  generalized  Kirchhoff- Love  shell  theory  to  thin  cylindrical  shells.  A  more  detailed 
shell  theory  uses  Reissner  (3)  and  Mindlin  (4)  (RM)  assumptions,  which  relaxes  the 
restriction  that  plane  sections  remain  normal  to  the  midplane.  A  nonlinear  shell 
theory,  the  Simplified  Large  Rotation  (SLR)  theory,  was  developed  by  Palazotto 
and  Dennis  (5) .  This  theory  has  been  validated  and  used  for  a  variety  of  loading 
conditions,  and  is  the  theory  used  in  this  research.  The  power  of  this  theory  is  the 
ability  to  account  for  the  geometric  nonlinearities  that  occur,  especially  for  a  thick 
shell  or  for  instability  such  as  buckling.  It  has  improvements  over  both  the  Donnell 
and  RM  theories,  as  will  be  described  in  Chapter  2. 
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Hatfield  (6)  used  the  SLR  theory  and  Finite  Element  Methods  (FEM)  to  ana¬ 
lyze  composite  cylindrical  shells  with  a  large,  rectangular  cutout  under  compression 
loading.  Various  configurations  were  chosen  and  compared  against  experimental  re¬ 
sults.  Both  quasi-isotropic  and  cross-ply  lay-ups  were  studied,  in  thicknesses  ranging 
from  8  plies  [1.04  mm  (0.041  in)]  to  24  plies  [3.13  mm  (0.123  in)].  Quasi-isotropic 
refers  to  combinations  of  0°,90°,  and  ±45°  laminates,  to  account  for  longitudinal, 
transverse,  and  shear  loads  (hence,  responding  similar  to  an  isotropic  material). 
Cross-ply  refers  to  a  layup  containing  0°  and  90°  oriented  plies.  This  work  studied 
the  geometric  nonlinearities  present  during  large  displacements  and  moderately  large 
rotations — such  as  during  collapse — as  well  as  the  effects  of  transverse  shear  in  thick 
shells. 

1.2.2  Cutouts.  Senocak  and  Waas  (7)  pursued  optimum  reinforce¬ 
ment  for  a  laminated  cylindrical  shell  with  a  cutout.  Collapse  due  to  a  static  buckling 
load  was  determined  using  the  Donnell  shell  theory.  Then  a  stress  function  was  as¬ 
sumed  and  applied  to  the  Donnell  equations.  This  resulted  in  an  expression  for  the 
forces,  bending,  and  twisting  moments  if  a  cutout  were  present.  Once  this  expres¬ 
sion  is  known,  the  forces  and  moments  required  to  resist  bending  can  be  determined. 
From  this  information,  a  reinforcement  is  designed  to  increase  the  shell’s  collapse 
load. 

Hilburger,  Waas,  and  Starnes  (8)  studied  the  interaction  between  bending  and 
membrane  reactions  due  to  the  presence  of  a  cutout  in  a  cylindrical  shell  under  a 
static  compression  and  internal  pressure  load.  The  nonlinear  code  STAGS  was  used 
in  the  analysis,  and  it  was  determined  the  cutout  introduces  local  buckling  which, 
reduces  the  overall  collapse  load  of  the  cylinder. 

Kapania,  Haryadi,  and  Haftka  (9)  used  an  alternating  global/local  analytical 
approach  to  determine  the  stress  field  around  circular  and  elliptical  holes  in  flat  com¬ 
posite  plates.  First,  a  global,  numerical  solution  is  obtained  using  the  Ritz  method 
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with  an  assumed  displacement  function.  These  resulting  displacement  and  slope 
boundary  conditions  are  then  applied  to  a  local  finite  element  model  to  determine 
the  stress  field  around  the  cutout.  For  a  circular  hole,  the  local  analysis  proved  to  be 
extraneous,  while  for  elliptical  holes,  the  extra  local  analysis  increased  the  accuracy 
of  results. 

Rhee,  He,  and  Rowlands  (10)  take  a  similar  approach  to  determining  the  stress 
around  a  cutout  in  a  composite.  First,  a  finite  element  analysis  is  performed  on 
a  plate  without  a  cutout,  where  a  single  element  encompasses  the  area  normally 
occupied  by  the  cutout.  The  displacement  boundary  conditions  for  this  element  are 
then  applied  to  a  Gerhardt  hybrid  element  (11).  This  hybrid  element  contains  a 
cutout,  which  must  be  described  by  a  mapping  function.  The  displacement  around 
the  cutout  are  determined  from  the  boundary  conditions  and  the  cutout  geometry, 
which  can  be  related  to  photomechanical  fields  such  as  Moire’s.  This  information 
can  then  be  used  to  determine  the  stress. 

Numerical  stress  distributions  for  a  thick  (nonlinear)  curved  shell  with  circular 
or  elliptical  cutouts  was  developed  by  Chao,  Diankui,  Xingrui,  and  Benli  (12).  The 
complex  variable  method  and  conformal  mapping  are  used  to  study  stress  concen¬ 
trations  in  thick  shells.  The  complex  stress  function  converts  the  shell’s  governing 
equations  into  three  Helmholtz’s  equations.  The  stress  function  is  represented  as  an 
infinite  series. 

1.2.3  Failure  Criteria.  Another  approach  to  analyzing  com¬ 
posite  materials  is  through  the  use  of  failure  criteria.  In  general,  a  failure  criterion 
compares  load  parameters,  such  as  the  current  state  of  stress,  with  the  material’s 
determined  parameters,  such  as  ultimate  strength  values.  If  the  criterion  exceeds 
this  value,  failure  is  determined  to  occur.  One  of  the  early  developers  of  a  failure  cri¬ 
terion  for  composites  was  Tsai  (13).  His  work  was  incorporated  with  other  theories 
to  develop  such  criteria  as  Tsai-Hill  (14)  and  Tsai-Wu  (15).  The  problem  with  the 
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Tsai-Hill  criterion  is  its  inability  to  account  for  differing  failure  strengths  in  tension 
and  compression  (16).  Hoffman  (17)  corrected  this  problem  by  adding  a  quadratic 
expression.  Tsai-Wu  used  a  general  quadratic  expression  to  add  robustness,  but 
this  resulted  in  a  cross-ply  material  characteristic  which  is  not  easily  determined 
(16).  Hashin  (16)  developed  a  simple  but  powerful  criterion  tailored  to  composite 
materials.  It  can  handle  differing  strengths,  and  determines  whether  failure  is  due 
to  tension,  compression,  or  shear  in  either  the  fiber,  matrix,  or  by  delamination  be¬ 
tween  layers.  This  is  the  criteria  used  for  the  present  research  and  will  be  discussed 
in  more  detail  later. 

Lee  (18)  incorporates  failure  criteria  with  finite  element  analysis  to  1)  deter¬ 
mine  the  stress  distribution  2)  identify  damage  and  failure  modes  and  3)  compile 
damage  accumulation  to  determine  ultimate  strength.  The  failure  criterion  devel¬ 
oped  by  Lee  is  similar  to  Hashin’s,  however,  it  does  not  incorporate  as  many  terms. 
Lee  identifies  procedures  for  a  progressive  failure  approach  in  a  finite  element  anal¬ 
ysis.  First,  the  finite  element  analysis  is  performed  and  the  stress  distribution  de¬ 
termined.  Once  the  stress  is  known,  the  failure  criterion  is  applied.  If  damage  has 
occurred,  the  appropriate  stiffness  terms  are  reduced.  The  externally  applied  stress 
is  then  increased  and  the  process  repeated.  This  continues  until  failure  of  the  entire 
panel  occurs.  This  approach  is  used  as  the  basis  in  this  research.  It  is  important  to 
note  that  Lee  admits  no  delamination  occurred  during  his  analysis,  which  does  not 
seem  to  match  experimental  results. 

Falzon,  Steven,  and  Xie  (19)  used  the  Tsai-Hill  failure  criterion  to  optimize 
cutouts  in  composites.  They  performed  a  finite  element  analysis  where  elements 
were  eliminated  in  the  regions  in  which  the  highest  failure  occurred.  This  allowed  a 
broader  redistribution  of  stress,  until  the  cutout  shape  was  optimized.  As  expected, 
for  an  equal  cross-ply  load,  the  optimal  cutout  is  a  circle.  For  a  pure  shear  load,  the 
optimal  cutout  is  square  rotated  45  degrees  from  the  panel  edges.  Under  combined 
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cross-ply  and  shear  loads,  an  elliptical  shape  is  optimal.  The  basis  for  the  present 
research  explores  geometries  which  are  not  allowed  to  be  optimal. 

Spottswood  (20)  used  failure  criterion  in  a  progressive  damage  model  to  de¬ 
termine  failure  in  a  cylindrical  shell  under  a  transverse  load.  Three  failure  criteria, 
the  maximum  stress  criterion,  Lee’s  criterion,  and  Hashin’s  criterion  were  used  in 
comparison. 

1.2.4  Fatigue.  A  fatigue  load  is  a  load  where  the  applied  magnitude 
is  below  the  failure  threshold  but  is  applied  over  many  cycles.  This  cyclic  load  causes 
micromechanical  damage  to  the  structure  (21).  This  minute  damage  causes  local 
stress  concentrations,  which  locally  increase  the  damage  zone.  Finally,  over  time, 
this  micro-damage  coalesces  into  larger,  macro-damage.  Once  the  macro-damage 
is  present,  stress  concentrations  occur  over  much  larger  regions.  This  causes  the 
macro-damage  to  grow  and  expand  to  an  increasingly  larger  region,  until  finally,  the 
structure  can  no  longer  carry  the  static  load. 

The  configuration  or  geometry  of  a  structure  also  affects  its  fatigue  life.  Both 
research  and  experience  have  shown  that  notches,  holes,  and  cutouts  cause  a  local 
increase  in  the  stress  distribution  (21).  Research  has  also  shown  that  sharper  notches 
and  corners  cause  a  larger  local  stress  concentration  than  more  rounded,  smooth 
cutouts  (21).  This  is  comparable  to  thinking  of  stress  as  water  flow  and  the  structure 
as  an  aquaduct.  A  sharp  reduction  in  the  aquaduct  size  requires  the  water  to  increase 
its  velocity  almost  instantaneously,  resulting  in  a  large  disruption  of  water  flow  in  the 
region  of  reduction,  while  a  gradual,  smooth  reduction  allows  the  water  to  accelerate 
smoothly,  over  a  longer  distance.  Likewise,  sharp  cutouts  cause  a  disruption  of  stress, 
resulting  in  a  large  stress  concentration.  Obviously,  then,  it  is  favorable  to  eliminate 
cutouts  or  design  them  as  non-interfering  as  possible,  such  as  the  work  by  Falzon, 
et  al.  mentioned  above.  This  is  often  difficult  in  applied  structures,  however,  for 
reasons  such  as  geometry  limitations  or  manufacturing  difficulty. 
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Tools  have  been  developed  to  predict  the  time  until  a  material  fails  in  fatigue, 
or  its  fatigue  life  (21).  On  the  micro-level,  the  field  of  fracture  mechanics  attempts 
to  define  the  stress  concentrations  around  geometry  such  as  notches,  holes,  and  crack 
tips.  On  the  macro-level,  most  isotropic  materials  and  several  composite  laminate 
materials  have  fatigue  life  or  S-N  curves.  The  fatigue  stress  is  usually  denoted  S 
while  the  number  of  cycles  to  failure  is  N.  However,  most  damage  initiates  on  the 
micro-level  at  some  inherent  defect  which  is  impossible  to  predict.  Also,  the  number 
of  other  inherent  defects  within  the  material  affect  how  fast  the  damage  changes  to 
the  macro-level.  Therefore,  predicting  fatigue  life — either  through  experimental  or 
computational  means — is  not  an  exact  science. 

1.2.5  Fatigue  in  Composites.  Fracture  mechanics  and  other 
fatigue  analysis  techniques  were  initially  developed  for  isotropic  materials  (21).  For 
an  isotropic  material,  loading  results  in  a  single  crack  that  grows  until  the  structure 
can  no  longer  carry  the  applied  load.  A  composite  laminae,  made  up  of  fibers  and  a 
matrix,  are  stacked  into  a  laminate.  There  are  many  failure  modes  for  a  composite 
material,  and  they  may  develop  independently  or  as  a  result  of  other  types  of  failure 
(22).  There  are  also  several  mechanisms  to  resist  failure  growth.  If  the  matrix 
cracks,  crack  growth  is  resisted  by  the  continuous  fiber.  If  a  fiber  fails,  the  matrix 
can  distribute  the  failed  fiber’s  load  to  other  fibers.  If  a  ply  disbonds,  each  layer  is 
still  able  to  carry  in-plane  loads.  Therefore,  the  initiation  and  interaction  of  these 
damage  modes  must  be  considered  for  composite  laminate  fatigue. 

Agarwal  and  Broutman  (22)  identify  the  various  levels  of  fatigue  damage  as 
it  occurs  in  multidirectional  composites.  The  first  effect  is  the  debonding  of  the 
fiber  and  matrix  in  the  most  transverse  plies.  For  example,  the  plies  oriented  90° 
to  the  load  (and  hence  the  most  transverse)  would  develop  these  disbonds  first. 
This  develops  into  a  crack  which  extends  vertically  to  the  top  and  bottom  of  the 
ply.  These  cracks  arrest  at  plies  with  fiber  orientations  parallel  to  the  load  (or  the 
principal  direction).  These  cracks  occur  at  several  locations  within  the  transverse 
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plies,  until  they  saturate.  This  has  been  termed  the  characteristic  damage  state  or 
CDS  by  Reifsnider,  Henneke,  and  Stinchcomb  (23).  The  crack  tips,  however,  cause  a 
local  stress  concentration  at  the  interface  between  plies,  which  cause  delamination. 
This  delamination  grows  until  large  stresses  develop  in  the  principal  direction  plies, 
causing  matrix  cracking  and  fiber  debonding.  Also,  the  plies  are  separated  and 
become  independent,  further  increasing  local  stress  until  fiber  breakage  occurs.  Large 
scale  fiber- breakage  occurs  near  90%  of  the  composite’s  fatigue  life  (22) . 

Fatigue  of  composite  materials  was  also  studied  by  Talreja  (24).  He  proposed 
that  fatigue  of  composites  be  characterized  by  strain  instead  of  stress,  since  the 
state  of  stress  within  a  composite  can  vary  depending  on  ply  layup.  He  identifies 
three  important  regions  in  the  fatigue  life  of  a  unidirectional  composite.  The  first 
is  the  endurance  limit,  which  may  be  defined  as  the  value  of  strain  below  which 
fatigue  failure  does  not  occur.  Typically,  this  is  the  strain  below  which  matrix 
cracking  does  not  occur,  denoted  as  em.  The  upper  region  is  the  strain  where  fiber 
breakage  occurs,  or  ec.  For  all  practical  purposes,  this  is  the  same  as  static  failure, 
and  fatigue  life  is  minimal.  In  the  middle  is  the  region  of  fatigue,  where  matrix 
cracks  initiate.  This  damage  leads  to  further  cracks  and  delaminations,  which  are 
the  damage  growth  mechanisms  leading  to  fatigue  failure.  For  a  multidirectional 
laminate,  a  different  lower  bound  may  occur.  For  laminae  in  directions  transverse 
to  the  load,  disbonding  between  the  fiber  and  the  matrix  may  be  present.  This  is 
extreme  at  the  90°  orientation,  and  lessens  for  decreasing  fiber  orientation  angles. 
This  research  correlates  well  with  the  failure  progression  described  by  Agarwal  and 
Broutman,  above. 

One  of  the  fracture  mechanics  techniques  applied  to  composites  is  through 
determining  the  strain  energy  release  rate,  G.  This  method  has  been  used  by  Ericson, 
et.  al.  (25)  as  well  as  Jen  and  Sun  (26).  Another  method  is  a  strength  based 
approach.  In  this  approach,  an  expression  for  the  residual  strength  of  the  composite 
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is  developed  as  a  function  of  the  number  of  cycles.  According  to  Sendeckyj,  three 
key  assumptions  determine  this  approach  (27): 

1.  A  Weibull  distribution  describes  the  static  strength. 

2.  The  residual  strength  after  n  cycles  is  related  to  the  initial  static  strength 
through  a  deterministic  equation,  based  on  the  S-N  curve. 

3.  Failure  occurs  when  the  residual  strength— the  remaining  material  strength  of 
the  structure-reduces  to  the  value  of  the  applied  stress. 

This  method  has  been  explored  more  recently  by  Reifsnider  (28),  as  well  as  Schaff 
and  Davidson  (29). 

Stinchcomb  and  Bakis  (30),  reflecting  on  the  work  of  Jamison  and  Reifsnider 
(31),  note  that  the  stiffness  of  a  graphite  epoxy  laminate  reduces  over  the  fatigue 
life  of  a  composite.  An  equation  can  be  used  to  account  for  this  stiffness  reduction 
over  time.  This  technique  is  explained  by  Sendeckyj  (27).  Poursartip,  Ashby,  and 
Beaumont  (32)  used  this  technique  to  study  delamination  under  fatigue  loading,  and 
their  research  is  based  on  the  work  of  O’Brien  (33).  O’Brien  used  a  rule  of  mixtures 
to  account  for  the  loss  of  stiffness  due  to  delamination. 

Jen,  Kao,  and  Hsu  (34)  also  studied  the  onset  of  delamination  by  analyzing 
stiffness  reduction.  They  used  a  mathematical  model  for  the  stress  distribution,  and 
applied  the  resulting  stress  to  a  failure  criteria.  Once  the  criteria  was  exceeded,  it  was 
determined  that  delamination  occurred.  They  were  unable,  however,  to  determine 
the  resulting  stress  redistribution  due  to  the  delamination. 

Barboni,  Carbonaro,  and  Gaudenzi  (35)  also  studied  the  onset  of  delamination 
due  to  fatigue.  Using  a  Multi-layer  Higher-Order  Laminate  Theory  which  takes  a  two 
dimensional  approach  to  a  three  dimensional  problem,  they  characterize  the  out-of¬ 
plane  stress  which  develops  near  a  free  edge.  An  interlaminar  failure  criterion  is  used 
to  determine  when  and  where  the  onset  of  delamination  occurs.  Their  theoretical 
work,  however,  only  focused  on  the  onset  of  delamination.  Experimentation  was 
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also  performed,  and  they  concluded  delamination  is  not  critical  for  [0/90]  layups,  in 
static  or  fatigue  loading.  For  [+45/  —  45]  layups,  delamination  initiates  in  the  static 
case  near  the  panel  failure  value.  But  for  fatigue,  delamination  can  be  a  driver  in 
the  reduction  of  fatigue  life  for  this  layup. 

Another  approach  is  to  incorporate  fatigue  parameters  into  the  failure  criteria. 
This  approach  has  been  taken  by  Hashin  and  Rotem  (36),  Fawaz  and  Ellyin  (37), 
and  most  recently,  Philippidis  and  Vassilopoulos  (38).  These  criteria  are  similar  to 
quasi-static  failure  criteria,  except  they  contain  factors  that  are  functions  of  fatigue 
data.  This  data  is  normally  acquired  through  the  determination  of  an  S-N  curve. 
The  downside  to  this  approach  is  the  reliance  on  experimental  data.  The  amount  of 
required  data  can  be  extensive  for  different  configurations  and  loading  conditions. 

Work  by  Tahiri,  Henaff-Gardin,  and  Lafarie-Frenot  (39)  shows  some  interesting 
relationships  between  static  and  fatigue  damage  development  and  stiffness  reduction. 
Their  research  involved  T300/914  graphite  epoxy,  [+452/ —  452]2s>  specimens.  Under 
cyclic  tests,  they  show  a  steady  increase  in  the  number  of  cracks,  while  for  quasi¬ 
static  loading,  cracks  do  not  develop  until  the  load  is  very  high.  Stiffness  reduction, 
however,  has  the  inverse  relationship.  The  material  stiffness  degrades  a  large  value  in 
quasi-static  loading  (up  to  75%),  but  only  varies  slightly  in  fatigue  loading  (around 
10%). 

As  seen  from  this  literary  review,  cutouts  and  fatigue  in  composites  are  all  being 
actively  studied,  with  many  different  approaches  being  taken.  Further  research  is 
warranted  in  these  areas  to  better  define  the  affects  on  composite  structures. 


1.3  Objective 

The  objective  of  this  research  is  to  develop  a  technique  to  characterize  the  dis¬ 
placement,  stress,  and  subsequent  failure  that  occurs  in  a  cylindrical  composite  shell 
with  a  square  cutout  under  a  fatigue  load.  This  research  involves  a  complex  combi- 
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nation  of  geometry,  materials,  and  loading.  Due  to  recent  advances  in  both  theory 
and  computing  power,  it  is  now  possible  to  analyze  such  complicated  combinations. 


1.4  Approach 

The  approach  of  this  research  is  to  use  a  Finite  Element  Model  (FEM)  to 
determine  the  stress  distribution  in  a  cylindrical  shell  with  a  large  rectangular  cutout. 
First,  a  static  analysis  is  performed  to  determine  the  shell  ultimate  load.  Then,  a 
fatigue  load  is  applied  at  a  percentage  of  the  ultimate  load.  This  is  done  for  tension 
as  well  as  compression. 

The  finite  element  code  was  developed  by  Dennis  and  Palazotto  and  uses  the 
SLR  theory  (5).  FEM  methods  provide  the  stress,  strain,  and  displacement  for  vir¬ 
tually  any  location  at  any  ply  in  the  cylinder.  This  provides  a  detailed  description 
of  the  stress,  not  only  within  the  element,  but  on  a  ply-by-ply  basis  as  well.  This  is 
a  benefit  over  current  experimental  analysis  techniques  such  as  strain-gages,  photoe¬ 
lasticity,  Moire  diffraction,  and  Stress  Pattern  Analysis  by  measurement  of  Thermal 
Emission  (SPATE),  which  have  a  limited  ability  to  combine  stress  intensities  and 
distributions  at  any  location. 

Failure  criteria  and  stiffness  reduction  techniques  are  incorporated  to  account 
for  damage.  The  Hashin  failure  criterion  is  chosen  for  its  detailed  stress  expressions 
and  its  ability  to  distinguish  between  matrix,  fiber,  and  delamination  failures  in 
either  tension  or  compression.  Once  failure  is  determined  to  occur  by  the  failure 
criterion,  the  appropriate  elemental  stiffness  terms  are  reduced  according  to  the 
type  of  damage.  The  stiffness  terms  are  also  arbitrarily  reduced  each  loading  cycle 
to  represent  the  effects  of  fatigue  on  a  composite.  Therefore,  this  is  more  of  a  damage 
accumulation/stiffness  reduction  approach  to  fatigue  versus  the  traditional  fracture 
mechanics  approach. 


1-11 


II.  Theory 

It  is  important  to  discuss  the  theory  behind  the  techniques  used  in  this  research.  The 
basis  of  the  analysis  is  the  SLR  shell  theory,  which  allows  use  of  more  complicated 
geometries  than  just  a  flat  plate.  It  is  also  necessary  to  discuss  the  finite  element 
techniques  used  to  implement  the  theory.  The  basis  for  the  development  of  the 
Hashin  criteria  is  described.  Finally,  stiffness  reduction  is  presented  as  a  technique 
for  analyzing  composite  materials  undergoing  cyclic  loading. 


2.1  Simplified  Large  Rotation  Theory 

As  mentioned  before,  the  SLR  shell  theory  is  used  for  two  reasons:  1)  it  allows 
2-D  analysis  of  a  3-D  problem  2)  it  accurately  handles  geometric  nonlinearities  from 
the  response  of  a  curved  panel.  The  SLR  theory  is  fully  explained  in  (5).  For 
completeness,  many  of  the  important  aspects  of  the  theory  will  be  discussed  here. 

2.1.1  Assumptions.  One  of  the  main  differences  between  the 
numerous  shell  theories  is  the  assumptions  used.  The  main  assumptions  used  in 
the  SLR  theory  are  as  follows:  1)  The  shell  thickness  is  small  compared  to  the 
in-plane  dimensions.  2)  The  shell  response  is  characterized  based  on  its  original, 
undeformed  configuration.  This  is  known  as  a  “Lagrangian”  approach.  Therefore, 
the  internal  strain  energy  is  characterized  by  Green’s  strain  tensor  and  the  Second 
Piola-Kirchoff  stress  tensor.  3)  The  material  response  is  assumed  to  be  linear,  while 
geometric  nonlinearities  are  allowed  to  occur.  A  geometric  nonlinearity  results  from 
large  displacements  and  moderately  large  rotations.  4)  Because  the  shell  is  thin, 
stress  in  the  out-of-plane,  or  <t3,  direction  is  zero.  5)  The  through-the-thickness,  or 
transverse,  shearing  strains  are  zero  on  the  top  and  bottom  surfaces  but  are  parabolic 
between  the  two.  Therefore,  the  transverse  stresses  also  have  this  representation.  6) 
While  the  stresses  and  strains  are  developed  in  the  Lagrangian  coordinate  system, 
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Figure  2.1  Shell  Reference  Geometry  and  Coordinate  System 

the  constitutive  matrix  is  based  on  an  Eulerian  coordinate  system.  This  is  assumed 
to  be  an  accurate  relationship  if  the  strains  are  small,  even  with  large  displacements 
and  rotations.  7)  The  theory  assumes  that  an  individual  fiber-matrix  “cube”  is 
transversly  isotropic,  but  a  laminae  of  unidirectional  fibers  in  a  matrix  behaves 
orthotropically. 

2.1.2  Background.  There  are  several  key  elements  to  the  SLR 
theory  which  describe  its  uniqueness  and  usefulness.  Since  a  shell  has  curvature,  a 
curvilinear  coordinate  system  is  used.  The  coordinates  for  curvature  in  two  directions 
are  shown  in  Figure  2.1  (5).  For  this  work,  however,  a  cylinder  is  assumed.  Therefore, 
the  following  references  from  Figure  2.2  are  used  (5). 

=  x  ( Ri  =  oo) 


2-2 


Figure  2.2  Cylindrical  Shell  Reference  Geometry  and  Coordinate  System 

ft  =  «  (H2  =  B)  (2.1) 

C  =  2 

Green’s  strain  tensor  is  developed  and  reported  in  pp.  22-26  of  (5).  For  the  stress  and 
strain  tensors,  a  shorthand  notation  is  listed  in  Table  2.1.  A  modified  representation, 
based  on  the  shell  assumptions,  is  used  in  the  SLR  theory.  First,  it  is  assumed  £3 
is  approximately  zero.  This  assumption  is  valid  since  the  thickness  of  the  laminate 
will  not  vary  significantly  during  loading.  Because  of  the  thin  shell  assumptions,  it 
is  reasonable  to  assume  the  in-plane  stresses  and  strains  are  more  dominant  than  the 
transverse  stresses  and  strains.  Therefore,  the  full  Green’s  representation  of  the  in- 


STRESS 

Table  2.1  Shorthand  Tensor  Notation 

STRAIN  COORDINATES 

CTn  =  <Tl 

Oi  = 

X  =  1 

022  = 

^22  =  ^2 

S  =  2 

033  — 

^33  =  C3 

Z  =  3 

023  =  04 

2623  —  €4 

SZ  =  4 

0r13  =  05 

2£l3  =  ^5 

XZ  =  5 

0"12  =  06 

2ei2  =  e6 

XS  =  6 
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plane  strains,  (e1,€2,  and  e6),  are  developed,  while  the  linear  terms  for  the  transverse 
strains,  (£4  and  e5),  are  used.  This  abridged  representation  of  the  strains  yields  a 
breakdown  in  continuity,  especially  for  the  transverse  shear  strains.  However,  for 
problems  of  moderate  to  small  rotations,  the  continuity  equations  are  approximately 
met.  The  compatibility  equations  are  satisfied  for  the  linear  case  where  only  small 
displacement  and  rotations  exist  (5). 

For  the  SLR  theory,  as  for  the  Reissner-Mindlin  (RM)  (3), (4)  theories,  plane 
sections  are  not  required  to  remain  normal  to  the  midplane.  This  allows  both  a 
bending  rotation  and  a  transverse  shear  deflection.  The  out-of-plane  motion  of  the 
element  is  not  captured  merely  by  the  change  in  slope  of  the  w  component,  but  also 
by  an  additional  rotation  degree  of  freedom,  T.  This  rotation  exists  in  the  x  and  s 
directions. 

2.1.3  Kinematics.  As  mentioned,  the  transverse  strains  assume 
a  parabolic  relationship  through  the  thickness.  This  representation  accounts  for  an 
internal  transverse  shear  which  reduces  to  zero  on  the  upper  and  lower  surfaces. 
In  order  to  develop  these  strains  in  the  laminate,  it  is  important  to  capture  the 
proper  kinematics  through  the  displacement  terms.  Applying  the  SLR  kinematics 
to  a  cylindrical  shell  results  in  the  following  representation  for  the  u,  v,  and  w 
displacements  (5): 

4 

u(x,s,z)  =  u0  +  z^x  -  ^z3(^x  +  w,x) 

v{x,s,z)  =  V0(l  -  ^)  +  z^fs  -  -^z3(^s  +  w,s)  (2.2) 

w(x,s)  =  w 

where  h  is  the  total  laminate  thickness,  uq  and  Vo  are  measured  at  the  midplane,  and 
z  is  the  distance  from  the  midplane.  These  kinematics  were  derived  from  a  Taylor’s 
series  expansion  for  the  z-direction  about  the  midplane.  This  approach  was  first  used 
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for  a  shell  by  Bassett  (40).  The  differences  in  many  of  the  shell  and  plate  theories 
depends  on  the  chosen  truncation  of  this  infinite  series. 

Shear  locking  is  a  problem  in  RM  kinematics.  As  the  thickness  decreases, 
bending  dominates  over  the  shear.  In  RM  theories,  the  shear  and  bending  terms 
are  the  same  order,  which  results  in  a  disproportionate  representation.  The  shear 
terms  form  a  “penalty  matrix”  (41).  This  penalty  matrix  overconstrains  the  finite 
element  analysis,  causing  it  to  “lock.”  In  the  SLR  theory,  there  is  a  higher  order 
representation  of  the  shear  terms,  which  reduce  to  zero  as  the  thickness  decreases. 
Therefore,  the  bending  terms  are  allowed  to  dominate  for  a  thin  model  and  shear 
locking  is  avoided. 

2.1.4  Elasticity  relations.  The  above  kinematics  are  applied  to 
the  Green  strain-displacement  relations  listed  in  (5).  As  mentioned  above,  the  full 
Green  strain  representation  is  used  for  the  in-plane  strains,  the  transverse  strains  are 
linear  in  displacement,  and  e3  is  assumed  negligible.  The  in-plane  strains  become: 

ei  =  e®  +  ZK\  +  z2k\  +  zzk\  +  zak\  +  z^k\ 

€2  —  €2  T  ZK 2  +  z? k>2  T  zz  1^2  T  K2  T  (2-3) 

C6  =  +  ZKS  +  ^6  +  Z3  K6  +  Z4  +  Z&K\ 

where  (J=l,  2,  6)  are  the  respective  midplane  strains  and  the  k,j  (J=l,  2,  6), 
(1=1,  2,  3,  4,  6)  are  individual  strain  components  corresponding  to  the  Ith  power  of 
z  (yet,  the  Kj  are  not  raised  to  a  power).  These  values  are  shown  in  Appendix  A  of 
(5).  The  transverse  shearing  strains  are: 

e4  =  j-^-r(w>i+^,)(1-^-) 

R 

£S  =  (t», ,+*.)(! (2.4) 
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One  should  note  the  parabolic  through-the-thickness  representation  of  the  e4  and  e5 
strains.  Also,  the  transverse  strains  are  zero  at  the  top  and  bottom  surfaces,  (|,  ;y;). 

In  shorthand  notation,  the  strains  can  be  represented  by: 


e-  +  zpKip 


(2.5) 


where  i=l,  2,  6  and  p=l,  2,  3,  4,  5,  6,  7  (See  (5),  equation  (2.48-49)  for  derivation 
of  this  expression) . 

Finally,  all  the  necessary  degrees  of  freedom  have  been  introduced,  [ u ,  v,  w ,  w,x ,  w,s ,  $x,  ^s], 
where  the  comma-subscript  denotes  derivative  with  respect  to  the  subscript. 


The  constitutive  relationship  for  an  orthotropic  unidirectional  laminae  is  given 
in  (1)  as: 
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0 

0 

e4 

Cr5 

0 

0 

0 

0 

G55 

0 

£5 

0 

0 

0 

0 

0 

1 

CO 

C3 

^6 

/ 

(2.6) 


where  Cij  is  a  function  of  the  material  properties  E\,E2,  G\2,  G 13,  G23,  V12,  ^21, 
and  1/23 •  Ei  is  the  Young’s  modulus,  Gij  is  the  shear  modulus,  and  is  the  Poisson’s 
ratio,  all  in  the  respective  directions.  The  a3  component  is  assumed  to  be  zero.  This 
representation  allows  for  an  e3.  It  order  to  eliminate  this  term,  line  three  of  the 
above  equation  is  solved  for  e3: 


C 13  C23 

£3  —  “TT-6!  -  Tt- e2 
<-'33  <^33 


(2.7) 
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where 

A  =  1  -  ^12^21 

For  a  multidirectional  laminate,  the  constitutive  equations  for  each  ply  must 
be  transformed  into  the  global  coordinate  system.  This  is  done  through  the  following 
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transformation: 


/  \ 

Qn 

Ql2 

0 

0 

0 

/  > 

ei 

Ql2 

Q22 

0 

0 

0 

^2 

<76 

>  =  [T*] 

0 

0 

<2  66 

0 

0 

[T*]t  < 

^6  > 

cr4 

0 

0 

0 

Q  44 

0 

^4 

<75 

\  / 

0 

0 

0 

0 

Q55 

^5 

V.  / 

This  requires  a  tensoral  transformation.  Recall  from  Table  2.1,  the  current  shorthand 
representation  of  the  strain  is  not  tensoral.  Therefore,  the  strain  is  converted  into 
tensoral  strain,  transformed,  and  then  converted  back  to  shorthand  notation.  The 
transformation  matrix  [T*]  is  then  developed  by: 


pr)  =  [fl][T][/r] 


(2.11) 


where 


[R] 


1  0  0  0  0 
0  10  0  0 
0  0  2  0  0 
0  0  0  2  0 
0  0  0  0  2 


m  = 


c2  s 2  — 2cs  0  0 

s2  c2  2cs  0  0 

cs  —cs  (c2  —  s2)  0  0 

0  0  0  c2  -s2 

0  0  0  s2  c2 
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and 


[iT] 


1  0  0  0  0 

0  10  0  0 

0  0  I  0  0 

0  0  0  i  0 

0  0  0  0  I 


c  =  cos(<3>)  s  =  sin(<3>)  <l>=fiber  orientation 


The  equations  resulting  from  this  transformation  can  be  broken  into  in-plane 
and  transverse  portions: 


/  N 

$11  Q 12  $16 

/  \ 

el 

\  a2 

>  = 

$12  $22  $26 

< 

^2  > 

0% 

V  > 

$16  $26  $66 

^6 

V.  / 

and 


where  the  individual  terms  are: 


Q44  Q  45 
Q45  Q 55 


€4 

£5 


(2.12) 


(2.13) 


Qn  =  Qnc4  +  Q2254  +  2Q\2(?  s2,  +  4  Q^c2  s2 
Q22  ~  Qn$4  +  $22  c4  +  2Ql2c2s2  +  4  QqqC2s2 

Q12  =  Ql\C2 s2  +  Q22C? S2  +  Qwic4  S4)  ~~  (2-14) 

Qie  —  Qnc3s  —  Q22c^3  +  $i2(C53  —  c3s)  +  2Q^{cs^  —  c3  s) 

^26  “  QllC53  —  Q22(? $  +  Ql2(c3^  ”  C53)  +  2Qq$(c^S  —  CS3) 

$66  —  Q\\(?s2  +  Q22C2S2  —  2Q12cV  +  $66  (c2  —  S2)2 
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and 


Q  44  —  Q44C  "F  Q55S 

Q55  =  Q44s2  +  Q55c2  (2.15) 

Q45  Q44CS  Q55CS 

again, 

c  =  cos($)  s  =  sin($)  <3>=fiber  orientation 

The  transformed  constitutive  matrix  coefficients  are  denote  by  a  bar,  Q^.  Note 
that  the  transformations  are  on  a  ply-by-ply  basis,  since  each  ply  may  contain  a 
unique  orientation.  Later,  these  will  be  integrated  over  the  thickness. 

2.1.5  Strain  Energy.  The  potential  energy  of  a  shell  is  described 
as  the  internal  strain  energy  plus  the  product  of  the  external  force  and  the  respective 
displacements.  This  is  abbreviated  as: 

UP  =  U-V  (2.16) 

In  general,  V  can  be  determined  by 

V  =  PiUi 

where  Pi  is  a  force  vector  component  and  Ui  is  the  associated  displacement  compo¬ 
nent.  The  strain  energy  term  can  be  represented  by: 

V  =  \  jyO?  ■  <LjiV  (2.17) 

From  equation  (2.12),  cq  =  and  equation  (2.17)  can  be  written  in  terms  of 

strain  only.  For  convenience,  the  strain  energy  is  broken  into  in-plane  and  transverse 
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energies: 


U  =  U1  +  U2 

U\  =  -  [  [ {Qn((<l  + zp Kip)2 +  Q 22(^2  +  zPk^p)2 

z  JnJh 

2Ql2(el  +  ^Pk1p)(£2  +  ZT  K2r)  +  Q66(e6  +  ZP  Kep)2 

2^16  (el  +  2?>ft:lp)(e6  +  ZT  Ker) 

2Q26(4  +  zpK2p)(el  +  zrK6r))dzdQ  (2.18) 


U2  =  ~  [  [ {Qu(e4  +  Z2^)2  +  Q5?>(e5  +  Z2^)2 
z  Jn  Jh 

2Q45(e4  +  z2K42)(€^  +  z2  K^))dzdVt  (2.19) 

where  p,  r  =  1,  2,  3,  4,  5,  6,  7. 

The  evaluation  of  this  expression  is  lengthy,  and  only  the  in-plane  equations 
are  shown.  Similar  terms  are  grouped  as  follows: 

U\  =  —  [  {u\  +  U2  +  u$)dVl  (2.20) 

2  Jn 

where 


(2.21) 

(2.22) 
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+  2  (ft jiftj4  +  ftj 2 ^3 ) ij  ~t”  (2ftj]_ftj5  "f”  2ftj2^i4  d~  ^z3 ) 

+  2(ftjifti6  +  Kj2Ki5  +  ftj3^z4)A  +  (2/^jl^z7 

+  2ftj2^z6  +  2  ft  j  3/^5  +  ^j4^i4:)Jij  (2.23) 

+  2(ftj2^i"7  +  ftj3^z6  +  ^74^5)  A7  +  (2ftj3^z7 
+  2 ft j 4 ft i 6  ~t~  K'jS^'ih ) -^z?  ~f”  ( 2 ftjf 4 ft z 7  "t~  5 ^z6 ) -^ij 

+  2ft^5ftj7  “h  ftj‘6^z6 ) B{j  d”  2 K>j 6 ft z 7 ^zj  d~ 


and  i,j  =  1,  2,  6  and  p,  r  =  1, 2, 3, 4,  5,  6,  7. 

The  elasticity  arrays  (A^Aj,  etc.)  are  determined  when  the  constitutive 
matrices  are  multiplied  by  different  powers  of  2:  and  integrated  over  the  thickness, 
yielding: 


[Aj  5  Aj ,  Dij ,  Eij ,  Fij ,  Gij ,  Hij ,  I{j ,  J{j ,  Kij ,  Aj ,  A? ,  -Aj  ,  Aj ,  A?] 

[  Qy  [1,  z,  z2,  z3,  z4,  z5,  z6,  z7,  z8,  z9,  z10,  z11,  z12,  z13,  z14]dz  (2.24) 
Jh 


where  i,  j  =  1,2,  and  6.  The  first  three  terms  are  identical  to  the  Aj,  Aj,  and  Aj 
matrices  from  classical  laminate  theory  (42).  The  same  procedures  for  the  transverse 
equations  produce: 


U2  * —  ~  f  (^ra^n^mn  H"  26nftm2-Dmrz  T  ^n2^m2  An)^  (2.25) 

*  </n 

and 

mn,  Ann,  ftm]  =  [  Qmn[l,Z2,  ZA]dz  (2.26) 

Jh 

where  m,  n  =  4,  5.  The  elasticity  arrays  comprised  of  odd  powers  of  2  (e.g.  By,  E^, 
,  . . . ,  Sij)  are  eliminated  if  a  symmetrical  ply  layup  is  assumed. 

In  FEM,  the  above  integrations  are  performed  numerically  in  a  technique 
known  as  Gauss  quadrature  (41).  In  Gauss  quadrature,  the  integration  is  evalu- 
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ated  at  certain  points,  known  as  Gauss  points.  If  enough  Gauss  points  are  chosen, 
the  numerical  technique  will  accurately  approximate  the  actual  integration. 


2.2  Finite  Element  Method 

In  statics,  force  and  displacement  are  related  by  a  stiffness  coefficient.  The 
simplistic  equation  that  defines  this  relation  is: 


F  =  kd 


where  F  is  the  force,  k  is  the  stiffness,  and  d  is  the  displacement.  This  equation 
applies  to  a  spring  just  as  well  as  other  continuum.  In  the  previous  section,  the  con¬ 
stitutive  relationship  between  stress  and  strain  was  introduced.  Generally  speaking, 
stress  is  a  function  of  the  force,  and  strain  is  a  function  of  displacement.  Therefore, 
the  constitutive  matrix  is  analogous  to  the  stiffness  coefficient,  k. 

In  the  previous  section,  the  potential  energy  functional,  11^,  was  also  intro¬ 
duced.  It  was  comprised  of  the  strain  energy  (a  relation  of  stress  and  strain)  and 
the  applied  force  times  displacement.  Again,  the  force,  displacement,  and  stiffness 
terms  are  prevalent  throughout.  In  this  way,  the  potential  energy  of  the  continuum 
can  be  used  to  solve  the  response  to  displacements  and/or  forces.  This  problem  and 
others  can  be  represented  by  differential  equations  (known  as  the  strong  form)  or 
integral  equations  (known  as  the  weak  form)  (41).  The  conditions  of  the  strong  form 
must  be  met  everywhere,  while  conditions  of  the  weak  form  must  be  met  only  in  an 
average  sense. 

A  functional,  such  as  the  potential  energy,  is  an  integral  expression  of  the 
differential  equations  that  govern  a  problem.  Therefore,  a  functional  is  in  the  weak 
form.  For  a  functional  to  represent  a  problem,  it  must  satisfy  three  things: 

1.  Boundary  conditions. 
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2.  Compatibility. 

3.  Equilibrium. 

Examples  of  each  for  a  statics  problem  are  as  follows:  If  the  end  of  a  continuum 
is  clamped  to  a  rigid  structure,  the  functional  must  not  allow  answers  containing 
displacements  at  that  point.  This  is  an  example  of  a  boundary  condition.  Also,  a 
functional  that  yields  a  discontinuous  solution,  such  as  having  a  jump  or  separation, 
violates  compatibility  for  the  continuum.  If  equilibrium  is  not  met,  the  result  would 
imply  that  forces  are  not  balanced. 

For  a  simple  rod,  it  might  be  possible  to  directly  solve  the  equation  F  =  kd. 
Most  problems,  however,  are  just  too  complex.  Therefore,  the  geometry  of  the 
continuum  is  broken  into  discrete,  or  finite,  representations.  It  is  now  possible  to 
solve  for  each  discrete  element,  and  then  sum  the  response  over  the  entire  continuum. 
While  each  element  may  not  meet  all  the  necessary  criteria,  it  is  possible  to  meet 
the  above  conditions  on  an  average  sense  over  the  entire  continuum. 

2.2.1  Building  the  Finite  Element.  The  FEM  model  that 
represents  the  SLR  theory  is  a  36  degree-of-freedom  curved  rectangular  element, 
shown  in  Figure  2.3  (5).  The  corner,  or  vertex,  nodes  have  seven  degrees  of  freedom, 
[u,v,w,w,x  ,w,s  There  are  four  “mid-side”  nodes  with  degrees  of  freedom 

u  and  v.  In  order  to  define  the  degrees  of  freedom  along  the  edges  of  the  element,  the 
displacements  or  rotations  are  interpolated  between  the  nodes.  The  type  of  function 
used  for  interpolation  regulates  the  level  of  continuity  obtainable.  The  notation  Cm 
is  used,  where  the  mth  derivative  of  the  function  is  continuous  (41).  For  example, 
if  the  function  <j>  has  C°  continuity,  4>  is  continuous,  but  its  derivative,  ^  is  not. 
C°  continuity  is  sufficient  for  extension,  since  only  the  displacements  are  required 
to  be  continuous.  Bending  rotation  is  defined  by  a  slope,  which  is  developed  by 
a  derivative  of  the  w  degree-of-freedom.  Therefore,  C 1  continuity  is  required  for 
bending  deflection.  For  the  present  theory,  only  C°  continuity  is  required  for  the 
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displacement  degrees-of- freedom,  u,v,^x,  and  Degrees-of-freedom  w,w,x,w,s 
define  bending  deflection  and  require  C 1  continuity. 

The  degrees  of  freedom  are  interpolated  by  multiplying  a  matrix  of  shape  func¬ 
tions  by  an  array  containing  the  degrees  of  freedom.  The  shape  functions  must  meet 
the  boundary  conditions  at  the  nodes,  while  representing  the  proper  interpolation 
between  nodes.  There  are  various  types  of  shape  functions.  Since  only  C°  continu¬ 
ity  is  required  for  u,  v,  ^x,  \I>S,  a  linear  or  Lagrangian  shape  function  may  be  used. 
However,  since  the  u  and  v  displacements  for  a  shell  are  quite  complex,  it  is  best 
to  use  higher  order  shape  functions  (5).  A  quadratic  representation  provides  more 
accuracy.  For  w  and  its  slopes,  a  Hermitian  shape  function  imparts  C 1  continuity. 

A  finite  element  is  usually  established  in  its  own  internal  coordinate  system. 
The  relationship  between  the  finite  element  coordinate  system  and  the  global  system 
is  a  Jacobian  matrix,  J.  The  Jacobian  is  found  by  using  chain-rule  differentiation 
between  coordinate  systems  (41).  Now,  the  physical  displacements,  d,  can  be  written 
in  terms  of  the  elements  degrees  of  freedom,  q  (5): 

{■ d }  =  [r][D]{<,}  (2.27) 


where 


{d}T  =  [u,  U,  1  ,  U, 2  ,  V ,  V,1  ,  V,2  ,  w,  W,i  ,  W,2  ,  W,U  ,  W ,22  ,  W,12  ,^1,  ^1,1,  ^1,2,  ^2,  ^2,1,  ^2,2. 


FI  =  [J]'1 

[D]  =  Shape  functions 
=  36  DOF  element  array 

where  {<?}  is  a  36x1  array,  [D]  is  an  18x36  matrix,  [F]  is  an  18x18  matrix,  and  {d} 
is  an  18X1  array. 


2-16 


2.2.2  Applying  the  Potential  Energy  to  Finite  Element 
Analysis.  The  potential  energy  representation,  II,,,  for  the  SLR  theory  is  (5): 

np  =  9j{K  +  y  +  (2.28) 

where  K  is  a  matrix  of  constant  stiffness  coefficients,  Ni  is  a  matrix  of  stiffness  coef¬ 
ficients  that  are  a  function  of  the  displacement,  N2  is  a  matrix  of  stiffness  coefficients 
that  are  a  function  of  the  displacement  squared,  and  R  is  the  force  vector.  If  the 
stiffness  coefficients  are  divided  this  way,  the  matrices  are  easier  to  manipulate.  The 
first  variation  of  this  expression  with  respect  to  the  degree  of  freedom  array,  {<?}, 
using  the  calculus  of  variations,  yields: 

hllp  =  8qT[(K  +  —  +  ~^~)q  —  R]  =  0  (2.29) 

Since  8q  is  an  arbitrary  value,  and  the  trivial  case  5q  =  0  is  ignored,  the  equation  in 
the  bracket  must  equal  zero  to  satisfy  equilibrium: 

/V,  No 

{k  +  Y  +  T  ^~r  =  0  (2-30) 

Notice  this  is  the  simple  equation  for  statics  written  at  the  beginning  of  this  section, 
F  =  kd.  The  static  behavior  of  an  arbitrary  shell  has  now  been  represented  through 
finite  element  modeling. 

Solving  equation  (2.30)  is  not  easy,  since  N\  and  N2  vary  with  the  displacement, 
making  the  equation  nonlinear.  This  equation  is  solved  through  the  iterative  Newton- 
Raphson  (NR)  technique,  which  is  a  powerful  tool  in  computational  analysis  (41). 
For  shorthand,  let: 

N.  No 

(K+I-f  +  I-f)q-R  =  F(q) 
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A  small  increment  is  added  to  q  and  the  Taylor’s  series  expansion  of  F(q  +  A q)  =  0 
is  developed.  The  higher  order  terms  are  neglected  (since  A q  is  small),  yielding: 


If 


then 


dF 

F(q  +  A  q)  =  F(q)  +  -7^- Ag  +  ■  •  •  =  0 


dF 


N,  No 

F{q)  =  {KFJfF1-f)q-R 


dF 

—  =  (K  +  N1  +  N2) 


(2.31) 


(2.32) 


(2.33) 


remembering  that  Ni  =  f(q)  and  N2  =  f{q2)- 

(K  +  Ni  +  N2)  is  the  tangent  stiffness  matrix,  denoted  as  Kr ■  Equation  (2.32) 
can  now  be  written: 


Kx^q  =  — F(q ) 


(2.34) 


The  current  values  of  q  are  substituted  into  N\  and  N2  of  Kt ,  and  equation  (2.34) 
is  solved  for  A q.  This  value  is  added  to  q  and  the  right-hand  side  (RHS)  is  updated. 
Remember  that 


Ni  No 

F(q)  —  {K  +  ~~  H  ~ )q  ~  R 


For  the  exact  solution,  F(q)  would  equal  zero.  For  FEM,  the  process  is  repeated 
until  F(q)  reduces  to  a  small  value.  This  is  done  through  a  convergence  check.  While 
there  are  several  ways  to  calculate  convergence,  the  chosen  method  for  the  SLR  FEM 
is: 


VTMF-VTM-ir' 


x  100%  <  TOL 


(2.35) 


VEMT 

where  i  is  the  degree  of  freedom  number,  r  is  the  increment,  and  TOL  is  the  tolerance 
value.  A  typical  tolerance  value  is  0.1%.  Once  the  tolerance  is  reached,  the  next 
increment  of  load  or  displacement  is  applied. 
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At  this  point  it  is  important  to  note  the  difference  between  the  increment, 
A q,  and  the  iteration.  The  problem  is  solved  through  incremental  additions  of  the 
degrees  of  freedom,  q.  For  each  increment,  equation  (2.34)  is  iterated  until  tolerance 
is  met.  Once  tolerance  is  met,  the  next  increment  is  applied.  In  the  Newton-Raphson 
approach,  the  tangent  stiffness,  KT,  (which  is  the  slope  of  equilibrium  curve  at 
the  point  in  question)  is  used  to  project  forward.  Other  approaches  use  the  initial 
stiffness  or  the  secant  stiffness.  The  initial  stiffness  is  the  linear  stiffness  and  the 
secant  stiffness  is  the  slope  of  the  line  from  the  origin  to  the  point  in  question.  These 
approaches  do  not  require  calculation  of  the  tangent  stiffness  each  iteration,  which 
is  computationally  more  efficient.  However,  the  initial  and  the  secant  stiffnesses  do 
not  take  as  large  a  step  each  iteration,  requiring  more  iterations.  The  NR  approach 
reduces  the  number  of  iterations  required  for  convergence. 

It  is  also  important  to  note  the  described  procedure  is  displacement  control, 
where  the  displacements  (or  the  degrees  of  freedom,  q )  are  incremented.  It  is  also 
possible  to  increment  load,  R,  or  even  a  combination  of  the  two. 

2.3  Failure  Criteria 

Failure  criteria  are  a  comparison  of  applied  parameters  to  material  allowable 
parameters.  The  parameter  of  choice  is  usually  stress,  but  others  such  as  strain  or 
strain  energy  (43)  are  also  used.  According  to  Tsai  (44),  as  reported  by  Spottswood 
(20),  most  failure  criterion  are  classified  as  interactive  or  non-interactive.  A  non¬ 
interactive  criterion  does  not  consider  multiple  parameter  interactions,  examples  of 
which  are  the  maximum  stress  and  maximum  strain  theories.  A  non-inter  active 
criterion  is  best  suited  for  brittle  isotropic  materials  (20).  An  interactive  criterion 
attempts  to  account  for  load  combinations. 

The  ultimate  values  of  a  stress  failure  criterion  can  be  plotted  in  the  stress 
plane.  The  maximum,  or  constant,  stress  criterion  is  a  rectangular  block,  which 
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Figure  2.4  Failure  Surfaces  for  Various  Failure  Approximations 

overestimates  the  biaxial  stress  interactions  as  shown  in  Figure  2.4  (16).  A  linear 
criterion  will  underestimate  the  cross-ply  stress.  A  quadratic  curve  fit  is  a  better 
approximation  of  the  material  strengths,  and  higher  order  approximations  increase 
complexity  but  do  little  to  improve  the  fit.  A  smooth  quadratic  curve  works  well 
for  isotropic  materials,  but  difficulties  are  introduced  when  applied  to  composite 
materials  due  to  their  different  failure  modes  (16).  The  smooth  curve  does  not  dis¬ 
tinguish  between  different  failure  modes  and  will  therefore  result  in  an  incompatible 
expression  in  regions  of  failure  mode  interaction.  In  order  to  use  a  quadratic  failure 
expression,  and  yet  consider  the  different  failure  modes  of  composites,  a  piecewise 
smooth  criterion  was  developed  by  Hashin  (16).  Each  smooth  curve  represents  a 
composite  laminate  failure  mode. 

The  basis  for  the  Hashin  criteria  is  the  transverse  isotropic  properties  of  a  fiber- 
matrix  composite.  If  the  criterion  is  to  model  the  actual  material,  it  must  also  exhibit 
transverse  isotropy — that  is,  it  must  be  invariant  in  rotations  around  the  fiber,  or 
xi,  axis.  The  failure  criterion  is  thus  based  on  the  following  stress  invariants  for  a 
composite  material  (16): 

h  =  o-n 

/2  =  <722  +  CT  33 
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(2.36) 


h  —  a23  ~  a22(J33 

OT  Iz  =  “(<722  —  Vzij2  +  ^23 

r  ^.2  I  2 

h  —  al2+a13 

2  2 

h  =  2<7i20'23cr13  ~  <722<713  —  033<712 

In  the  transverse-isotropic  model,  the  1  direction  is  aligned  with  the  fiber.  Only 
quadratic  terms  are  desired,  and  the  cubic  term,  /5,  is  discarded.  The  resulting 
invariants  are  then  combined  to  form  the  general  representation  of  the  failure  criteria: 

A\Ii  +  B\Ii  +  A2I2  +  B2I2  +  C12I1I2  +  A^Iz  +  A4I4  =  1  (2.37) 


The  factors  A{,  Bi,  and  C ij  are  determined  depending  on  the  type  of  material  strength 
corresponding  to  the  invariants. 


2.3.1  Fiber  Failure.  Fiber  failure  in  tension  is  assumed  to  not 
depend  on  <722,033,  or  <123,  which  are  matrix  dominated  properties.  For  fiber  failure 
in  tension: 

1  1 
B\  —  — 2 —  and  — 

°FN  aFS 


where 


aFN  =  Ultimate  Fiber  Failure  in  Tension  Stress 
aFS  —  Ultimate  Fiber  Failure  in  Shear 


and  all  other  terms  are  assumed  to  be  zero.  This  leads  to  the  expression  where  fiber 
failure  in  tension  occurs  if: 


o_ n_ 
aFN 


+ 


(<7?2  +  Oil) 


>  1 


’ FS 


(2.38) 
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In  order  to  use  the  failure  criterion,  it  is  necessary  to  perform  experimental  tests  to 
determine  individual  failure  properties.  Fortunately,  these  values  have  been  deter¬ 
mined  for  most  composite  materials. 

2.3.2  Matrix  Failure.  It  is  assumed  that  the  fiber  carries  the 
longitudinal  stress,  an,  so  this  factor  is  not  present  in  the  matrix  failure  criterion. 
Equation  (2.37)  is  reduced  to: 


A2I2  +  B2I2  +  -d.3/3  +  A4I4  —  1 


(2.39) 


Matrix  failure  is  divided  into  two  types:  tension  and  compression.  For  the 
tension  case,  the  factors  in  equation  (2.39)  are  determined.  First,  A2  is  deleted  in 
favor  of  the  quadratic  term  B2.  Tensile  loading  for  the  matrix  occurs  in  either  the  2 
or  3  direction,  so 


B2  = 

As  = 

A4  = 


1 

aMNT 

1 

aMS 

1 

°FS 


and  the  matrix  tension  criterion  is: 


(022  +  033) 2  (023  —  022033)  .  (of 2  +  ^13) 


+ 


+ 


>  1 


a 


MNT 


'MS 


’ FS 


(2.40) 


The  compression  case  is  more  complicated.  First  a  simple  unidirectional  com¬ 
pression  loading  is  considered  and  an  expression  is  determined  in  terms  of  A2  and 
B2.  Then,  the  case  where  <r2 2  and  033  are  equally  applied  in  compression  is  analyzed 
(again,  ignoring  the  shearing  terms  of  J3  and  I4),  resulting  in  another  expression 
for  A2  and  B2.  These  two  variables  are  then  solved  for  by  these  two  equations. 
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In  the  resulting  expressions  for  these  two  variables,  the  terms  above  first-order  are 
discarded.  The  resulting  values  are  then  applied  to  equation  (2.39),  which  becomes: 


KSf)2  -  tK^B  +  gsa)  +  (022  -  ox)1  +  (<4  -  <^33)  +  (ft 2  +  dll  >  j  (2.41) 

&MNC  4 a\jS  aFS 

where 


oMnt  —  Ultimate  Matrix  Failure  in  Tension 
aMS  =  Ultimate  Matrix  Failure  in  Shear 
&mnc  —  Ultimate  Matrix  Failure  in  Compression 


2.3.3  Delamination.  Delamination  between  the  plies  can  also  be 
assessed.  The  criterion  is: 

(<4i  +  »a) 

°ls 

where 

a ds  =  Ultimate  Delamination  Shear  Strength 

This  does  not  include  the  transverse  normal  stress,  033,  which  should  be  incorporated 
for  a  more  accurate  determination  of  delamination.  Recall,  however,  that  in  the  SLR 
theory,  cr33  =  <j3,  which  is  assumed  to  be  zero.  This  factor  is  then  not  added  to  the 
failure  criteria. 

2.3.4  Maximum  Strain  Failure  Criterion.  The  need  for  a 
strain-based  failure  criterion  is  evident  from  the  fatigue  results,  discussed  in  Chapter 
4.  A  simple  strain-based  criterion  is  the  maximum  strain  failure  criterion.  For  this 
criterion,  failure  occurs  once  the  internal  strain  exceeds  the  material  ultimate  strain. 
For  the  fiber  direction: 

—  >  1  (2.43) 

6n 
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For  the  matrix  direction: 


(2.44) 


where 


1 


eN  =  Ultimate  Strain  in  the  1-Direction 
et  —  Ultimate  Strain  in  the  2-Direction 

While  not  as  exact  as  a  quadratic  criterion  such  as  Hashin,  it  does  identify  the  upper 
limit  of  strain  failure. 


2.4  Stiffness  Reduction 

Jamison  and  Reifsnider  (31)  identified  three  distinct  stages  of  stiffness  reduc¬ 
tion  in  a  graphite  epoxy,  [0/902]s,  laminate  as  shown  in  Figure  2.5.  Stinchcomb  and 
Bakis  (30)  comment  on  these  stages,  stating  that  Stage  I  is  an  rapid  initial  decrease 
primarily  due  to  matrix  cracking  and  some  minor  fiber  failure.  This  is  where  the 
characteristic  damage  state  (CDS)  occurs.  The  CDS  is  defined  as  matrix  cracking 
in  the  transverse  plies,  which  eventually  saturate  and  stabilize  (23).  Stage  II  is  de¬ 
scribed  by  Stinchcomb  and  Bakis  as  “an  intermediate  but  long  period  of  stiffness 
reduction  which  results  from  additional  matrix  cracking  in  off-axis  and  on-axis  plies, 
crack  coupling  along  ply  interfaces,  and  internal  delaminations”  (30) .  This  correlates 
with  the  fatigue  work  of  Talreja  (24),  also  noted  in  Chapter  1.  In  Figure  2.5,  it  is 
observed  that  during  Stage  II  the  stiffness  reduction  is  approximately  linear  with 
the  number  of  cycles.  Finally,  Stage  III  occurs  at  the  end  of  the  fatigue  life,  and 
is  due  to  large  delaminations  and/or  fiber  failures.  The  damage  mechanisms  are 
correlated  to  fatigue  life  through  Figure  2.6  (30).  Poursartip  et  al.  (45)  also  show 
similar  results  for  quasi-isotropic  materials,  especially  for  the  initial  Stage  I  drop  and 
the  approximately  linear  Stage  II  reduction  in  stiffness  with  number  of  cycles. 
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Figure  2.5  Stiffness  Reduction  Due  to  Fatigue  of  a  [0/902]s  Graphite  Epoxy  Lam¬ 
inate 


DAMAGE  MODES  DURING  FATIGUE  LIFE 


Figure  2.6  Damage  Development  During  Fatigue  Life  of  a  Composite  Laminate 
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Figure  2.7  Normalized  Stiffness  vs.  Fatigue  Life  for  High  and  Low  Fatigue  Loads 

Some  interesting  results  regarding  fatigue  and  stiffness  reduction  are  noted  by 
these  authors.  Their  experiments  involved  32-ply  T300/5208  graphite  epoxy  quasi¬ 
isotropic  laminates  (28).  For  fatigue  loads  near  45%  of  the  static  strength,  fatigue 
lives  were  on  the  order  of  450,000  cycles  to  runout  at  2,600,000  cycles.  Fatigue 
loads  near  55%  of  the  static  strength  resulted  in  fatigue  lives  only  around  100,000 
cycles.  Yet,  for  all  these  different  tests,  Stage  I  stiffness  reduction  occurs  through 
10-15%  of  the  fatigue  life  and  Stage  III  begins  at  90-95%  of  the  fatigue  life.  It  seems, 
then,  stiffness  reduction  is  a  favorable  parameter  to  use  in  characterizing  fatigue  in 
composites.  In  fact,  Stinchcomb  and  Bakis  (30)  show  in  Figure  2.7  that  at  both  low 
and  high  fatigue  loads,  the  change  in  the  stiffness  correlates  directly  to  percentage 
of  fatigue  life  (28). 

It  is  important  to  discuss  the  effects  loading  and  unloading  have  on  the  stress- 
strain  curve,  since  this  will  affect  the  stiffness.  Agarwal  and  Broutman  (22)  state 
if  a  composite  material  is  loaded  above  its  matrix  cracking  threshold,  a  change  in 
stiffness  will  occur  (22).  Unloading  occurs  along  a  straight  line,  but  one  with  a  slope 
different  than  the  initial  stiffness.  This  unloading  will  result  in  a  residual  strain, 
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much  like  an  elastic-plastic  model  for  an  isotropic  material.  Unlike  elastic-plastic 
modeling,  however,  this  residual  strain  will  reduce  to  near  zero  upon  increasing 
numbers  of  cycles.  This  result  is  also  proven  by  Tahiri  et  al.  (39)  in  plotting  <712  vs. 
€12-  They  show  a  build-up  of  residual  strain  for  quasi-static  loading  involving  two  or 
three  cycles.  Next,  the  initial  stress-strain  slope  is  compared  with  one  at  95%  of  the 
fatigue  life.  Both  slopes  originate  from  the  same  value  of  strain  (and  the  residual 
strain  is  no  longer  present),  but  the  latter  slope  is  reduced  from  the  initial  value  by 
about  19%. 
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III.  Numerical  Approach 

Finite  element  analysis  is  a  powerful  tool  to  solve  complicated  problems  by  breaking 
them  into  smaller,  easier-to-handle  elements.  There  are  several  techniques  which 
allow  each  analysis  to  be  tailored  to  the  problem  at  hand.  It  is  important,  however, 
that  the  finite  element  approach — at  least  on  a  global  or  average  sense — approximate 
the  physical  results.  For  this  research,  some  of  the  key  factors  in  the  finite  element 
analysis  are:  1)  the  boundary  conditions  and  inputs,  2)  application  of  the  failure 
criterion,  3)  iterating  on  damage,  4)  reducing  the  stiffness  to  account  for  fatigue, 
and  5)  the  finite  element  model. 


3.1  Boundary  Conditions  and  Inputs 

The  boundary  conditions  for  a  static  problem  are  its  constraints,  applied  dis¬ 
placements,  and  applied  tractions.  In  order  for  the  model  to  behave  as  the  actual 
structure  does,  the  proper  constraints  must  be  enforced.  Constraints  are  usually 
considered  degrees  of  freedom  that  are  fixed  at  zero.  In  the  Finite  Element  Method 
(FEM),  for  example,  if  an  edge  is  to  be  clamped,  then  the  degrees  of  freedom  of  the 
nodes  along  the  edge  are  forced  to  be  zero. 

The  inputs  are  what  cause  the  structure  to  displace,  and  hence  cause  internal 
strains  and  stresses  to  develop.  For  static  problems,  the  two  most  common  types  of 
inputs  are  displacements  and  forces.  The  forces  and  the  displacements  are  propor¬ 
tional  to  the  stresses  and  strains,  respectively.  Therefore,  an  equilibrium  relationship 
exists  between  force  and  displacement,  just  as  it  does  for  stress  and  strain.  For  linear 
problems,  once  one  is  applied  the  other  may  be  found  through  the  linear  relationship. 
As  mentioned  in  Chapter  2,  for  nonlinear  problems,  each  resulting  equilibrium  point 
is  based  on  the  previous  equilibrium  point,  and  an  incremental  approach  is  required. 
It  is  very  important  in  nonlinear  problems  whether  the  inputs  are  displacements, 
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Figure  3.1  Load  and  Displacement  Control 

forces,  or  a  combination  of  both,  as  the  equations  must  be  developed  in  a  way  to 
accommodate  the  type  of  input.  When  displacements  are  applied,  it  is  known  as 
displacement  control,  and  when  forces  are  applied,  it  is  known  as  load  control.  The 
necessity  to  use  one  approach  or  the  other  is  developed  when  two  values  may  result 
for  one  input.  This  is  evident  in  Figure  3.1.  For  an  applied  load  (the  load  control 
line),  there  are  two  equilibrium  points.  The  second  point  will  never  be  reached  by 
this  method.  The  displacement  control  line  easily  traverses  the  maximum  point  and 
accurately  determines  the  ultimate  load.  Combinations  of  load  and  displacement 
control,  such  as  the  Riks-Wempner  approach  (5),  are  much  more  complicated. 


3.2  Applying  the  Failure  Criterion 

The  Hashin  failure  criterion  is  easily  applied  to  the  finite  element  analysis,  as 
described  by  Spottswood  (20).  The  first  step  in  the  finite  element  analysis  is  the  ap¬ 
plication  of  the  boundary  conditions  and  inputs.  The  solution  is  then  iterated  using 
the  Newton- Raphson  technique  until  equilibrium  for  that  increment  is  reached.  The 
degrees  of  freedom,  (also  referred  to  as  the  displacements,  even  though  they  contain 
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rotations  as  well)  are  now  known  for  the  given  increment.  These  displacements  are 
converted  into  strain  using  the  strain-displacement  relations.  In  the  SLR  code,  this 
is  performed  in  the  STRESS  subroutine,  which  may  be  called  for  every  element.  The 
strain  is  determined  at  the  four  outermost  or  “corner”  Gauss  points  for  every  ply. 
This  is  done  through  nested  loops,  where  the  outer  loop  progresses  through  the  plies 
while  the  inner  loop  counts  through  the  four  Gauss  points  within  a  ply.  The  strains 
are  then  converted  into  stresses  through  the  constitutive  relationships,  again  inside 
the  nested  loops.  These  stresses  may  then  be  displayed  to  the  output  file  at  each 
Gauss  point  at  every  ply  in  every  element. 

The  material  strengths  for  the  given  composite  are  determined  in  the  local 
axis,  based  on  a  unidirectional  fiber  laminae.  The  stresses  must  be  transformed  into 
each  ply’s  local  axis  for  comparison  with  the  material  strengths  through  the  failure 
criterion.  A  permanent  flag  is  raised  if  the  failure  criterion  is  exceeded,  allowing  the 
accumulation  of  damage  through  the  increments. 

For  the  model  to  simulate  actual  damage,  it  must  become  more  flexible  as 
damage  occurs.  This  is  done  through  the  reduction  of  the  constitutive  terms.  Just 
as  terms  in  the  general  Hashin  criterion  are  tailored  to  the  component  evaluated,  the 
element  stiffness  terms  are  reduced  depending  on  the  type  of  failure.  Recall  equation 
(2.8): 
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(2.8) 


The  fiber  is  the  key  reinforcing  element  in  a  composite,  and  it  dominates  the 
properties  in  the  axial  direction.  In  the  material  axis  for  a  unidirectional  laminae, 
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the  1-direction  is  oriented  with  the  fiber.  Therefore,  fiber  failure  affects  all  proper¬ 
ties  that  involve  the  1  direction,  which  are  Qn, Q12, Q55,  and  Q 66-  These  terms  are 
reduced  if  fiber  failure  occurs.  This  is  for  tension  only,  since  fiber  failure  due  to  com¬ 
pression  is  a  complicated  failure  mechanism  and  dependent  on  the  micromechanical 
interaction  with  the  matrix  (16).  Fiber  failure  in  compression  is  not  considered. 

Matrix  failure  is  considered  in  tension  and  in  compression.  Once  failure  occurs, 
however,  the  result  is  the  same.  The  terms  Q22,Qi2,Qa,  and  C?66  are  reduced 
accordingly.  It  is  observed  that  failure  to  one  constituent  (e.g.  matrix)  does  not 
“zero  out”  an  element,  but  it  is  still  able  to  carry  load  through  the  other  constituent 
(e.g.  fiber).  This  is  one  benefit  of  composite  materials,  and  modeling  them  in  this 
way  is  more  appropriate. 

Application  of  the  delamination  failure  criterion  is  more  involved.  Delamina¬ 
tion  occurs  between  plies,  but  the  stresses  are  only  measured  in  the  middle  of  a  ply. 
Since  cr3  is  ignored  in  the  SLR  theory,  the  delamination  stresses  are  the  two  trans¬ 
verse  shear  stresses,  cr4  and  <75.  If  these  stresses  exceed  the  criterion,  delamination 
is  assumed  to  occur  between  the  plies.  This  is  an  approximation,  but  valid  within 
the  shell  theory’s  parabolic  representation  of  the  transverse  shearing  stress.  The 
parabolic  representation  ensures  that  if  the  stress  at  the  middle  of  the  ply  exceeds 
the  criterion,  the  stress  for  at  least  one  of  the  interfaces  will  also  exceed  the  criterion. 
This  is  true  for  a  symmetric  layup.  If  there  is  a  single  middle  ply  in  the  layup,  the 
parabolic  representation  is  a  maximum  at  the  middle  of  this  ply  and  this  reasoning 
breaks  down.  Even  this  case  is  considered  negligible  given  the  other  assumptions, 
such  as  <73  =  0.  As  the  shell  thickness  increases,  the  a3  component  becomes  more 
important.  A  code  that  accounts  for  the  03,  includes  the  delamination  effects  of  a 
free-edge,  and  determines  the  stresses  at  the  ply  interfaces  is  required  for  a  more 
accurate  determination  of  delamination  between  plies. 

Delamination  failure  is  modeled  by  reducing  the  Q44  and  Q55  stiffness  terms  in 
the  determining  ply.  If  Q44  and  Q55  are  eliminated  in  a  ply,  the  transverse  shearing 
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stress  cannot  transmit  across  this  ply.  This  is  similar  to  a  delamination  between 
plies,  which  also  prohibits  the  transverse  shearing  stress  from  crossing  between  plies. 


The  above  procedures  are  also  performed  in  the  STRESS  subroutine.  First,  the 
constitutive  terms,  Qij,  are  stored  in  dummy  variables,  AQ ^ ,  at  the  beginning  of  the 
subroutine.  Then,  in  the  inner  loop  of  the  STRESS  subroutine,  the  failure  criteria 
are  applied  at  each  Gauss  point.  If  a  criterion  is  exceeded,  the  appropriate  dummy 
term  is  reduced  to  a  small  number.  The  term  is  not  set  exactly  to  zero,  otherwise 
the  stiffness  matrix  would  become  singular  as  the  number  of  elements  with  damage 
increases.  The  arbitrary  factor  ^  is  used  in  the  program,  as  illustrated  in  this 
sample  reduction: 


AQij  = 


1 

4000 


AQij 


Each  AQi:j  term — whether  it  has  been  reduced  or  not — is  divided  by  four,  since  this 
process  occurs  for  each  of  the  four  Gauss  points  within  a  ply.  These  dummy  variables 
are  summed  into  another  dummy  variable  TQij,  which  is  initialized  to  zero  at  the 
beginning  of  the  STRESS  subroutine: 


TQ^  =  TQ^  + 


AQ^ 

4 


which  is  inside  the  inner  Gauss  point  loop.  If  no  failure  occurred,  the  rebuilt  con¬ 
stitutive  term  is  equivalent  to  its  original  value.  In  this  way,  the  constitutive  term 
is  reduced  proportional  to  the  amount  of  failure.  For  example,  if  failure  occurred  at 
two  of  the  four  Gauss  points  within  a  ply,  TQq  would  be  approximately  (ignoring 
the  residual  effect  from  the  multiplication)  50%  of  the  original  Qij  value.  This 
level  of  detail  provides  the  opportunity  for  an  accurate  representation  of  actual  com¬ 
posite  material  behavior.  As  before,  these  local  constitutive  terms  are  transformed 
into  global  coordinates  and  integrated  through  the  thickness. 
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3.3  Iterating  on  Damage 

In  nonlinear  finite  element  analysis,  an  incremental  approach  is  required.  For 
each  increment,  iterations  are  performed  until  equilibrium  is  reached.  This  incre¬ 
mental  approach  is  favorable  for  a  damage  accumulation  or  a  progressive  damage 
approach,  since  failure  is  evaluated  at  each  increment.  Failure  is  not  evaluated, 
however,  until  after  equilibrium  is  reached.  This  impact  is  shown  in  Figure  3.2. 
As  material  failure  occurs,  the  stiffness  of  the  material  decreases.  This  change  in 
stiffness  can  be  seen  in  the  figure  as  the  line  from  point  A  to  point  B.  The  FEM 
analysis  will  project  along  the  initial  stiffness  line  to  its  perceived  equilibrium  point 
at  C,  past  the  change  in  stiffness.  Then,  the  failure  criterion  is  implemented  and  the 
stiffness  terms  are  reduced. 

The  FEM  response  should  shift  to  the  actual  equilibrium  line.  This  shift  will 
vary  depending  on  the  type  of  input.  For  load  control,  the  same  increment  of  load 
will  result  in  a  larger  displacement,  or  line  CD  of  the  figure.  For  displacement  control, 
the  same  increment  of  displacement  will  result  in  a  lower  value  of  load,  which  is  line 
CE  of  Figure  3.2  .  This  is  known  as  “load  drop”  (20). 


Figure  3.2  Equilibrium  Responses  to  Failure  Depending  on  Type  of  Input 

Once  failure  occurs,  it  is  preferable  to  re-apply  the  same  increment  to  determine 
the  shift  to  the  new  equilibrium  value,  such  as  from  point  C  to  the  equilibrium  line 
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in  Figure  3.2.  If  not,  it  is  necessary  to  use  small  increments  to  reduce  the  amount  of 
error  at  each  determined  equilibrium  point.  This  requires  many  increments,  however, 
which  is  costly  in  terms  of  computational  time  and  negates  the  benefit  of  the  Newton- 
Raphson  approach.  Also,  this  approach  overestimates  the  final  failure  value  by  one 
step,  since  another  increment  is  required  to  produce  global  structure  failure  after  the 
“critical”  local  failure  occurs. 

Re-applying  the  same  increment  once  failure  occurs  is  termed  “iterating  on 
damage”.  This  requires  returning  to  the  previous  equilibrium  value  and  projecting 
forward  using  the  new  reduced  stiffness  matrices.  Previous  to  this  research,  this 
capability  did  not  exist  in  the  SLR  code.  A  flow  chart  of  the  algorithm  for  this 
process  is  shown  in  Figure  3.3.  Implementing  this  algorithm  requires  storing  the 
previous  equilibrium  position.  In  the  SLR  program,  this  is  done  through  storing  the 
global  displacement  values  at  the  previous  equilibrium  point.  As  the  failure  criterion 
is  applied,  a  flag  is  raised  if  new  failure  occurs.  This  flag  tells  the  code  to  re-apply 
the  same  increment,  starting  from  the  previous  equilibrium  value.  The  elements  are 
checked  for  failure  and  the  process  repeats  until  no  further  damage  occurs  at  that 
increment.  The  number  of  “iterations  of  damage”  depend  on  the  increment  size  and 
the  size  of  the  element.  The  program  then  progresses  to  the  next  increment.  In  this 
way,  the  increment  that  produces  global  structural  failure  will  be  identified.  Also, 
the  equilibrium  points  along  the  curve  will  be  accurate. 

The  effect  of  the  progressive  failure  on  an  element  is  shown  in  Figure  3.4,  where 
failure  in  an  individual  element  from  a  [0/90]2S  graphite/epoxy  panel  under  tension  is 
plotted  as  a  function  of  the  panel’s  normalized  displacement.  The  100%  value  is  the 
displacement  at  which  the  maximum  load  occurs.  An  element  near  the  cutout  was 
chosen  to  demonstrate  the  failure  progression  as  the  damage  area  grows  out  from  the 
cutout  corner.  The  failure  is  broken  into  matrix  and  fiber  failure  in  the  0°  and  90° 
plies  individually.  One-hundred  percent  failure  for  a  certain  ply  orientation  occurs 
when  all  the  Gauss  points  in  all  those  plies  exceed  the  failure  criterion.  It  is  noted 
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Figure  3.3  SLR  Algorithm  for  Iterating  on  Damage 
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Figure  3.4  Progression  of  Failure  in  a  Cross-Ply  Element  under  Tension 

that  fiber  failure  does  not  occur  in  the  90°  plies,  since  the  fibers  are  perpendicular 
to  the  primary  load.  Evidence  of  the  characteristic  damage  state  (CDS)  exists  in 
this  element  at  45%  of  the  ultimate  panel  value.  The  90°  plies  exhibit  69%  failure, 
which  then  stabilizes  until  90%  of  the  ultimate  value,  when  fibers  in  the  0°  direction 
fail.  Finally,  some  matrix  failure  does  occur  in  the  0°  plies,  but  it  is  obvious  from 
the  results  that  the  critical  constituent  in  the  0°  direction  is  the  fiber. 


3.4  Stiffness  Reduction  due  to  Fatigue 

The  fatigue  load  is  applied  at  a  fraction  of  the  static  failure  load.  If  the 
fatigue  value  is  above  the  matrix  cracking  threshold,  matrix  failure  as  well  as  possible 
localized  fiber  failure  will  occur.  The  new  equilibrium  value  is  found  and  the  load 
is  relaxed  to  its  initial  value,  completing  the  first  cycle.  The  shell  relaxes  to  a  value 
close  to  its  initial  value.  There  is  some  slight  residual  strain,  which  is  acceptable. 
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Without  further  adjustments,  the  shell  is  in  equilibrium.  No  further  failure  will  be 
produced  for  additional  cycles  between  these  two  load  values.  Fatigue  effects  must  be 
accounted  for  by  additional  steps.  This  is  done  by  reducing  the  material  properties, 
Ei,E2,  G\2 i  C?i3,  and  G2 3  for  the  entire  shell.  This  is  a  new  addition  to  the  SLR  code 
and  is  performed  in  the  subroutine  STRESS  at  the  beginning  of  each  cycle,  except  for 
the  first  cycle.  As  the  program  runs  through  the  stress  analysis  for  every  element,  the 
properties  Qn,  Qn,  Q22,  G12,  Gu,  and  G23  are  reduced  by  a  factor  (note:  Qn,Qn, 
and  Q22  are  all  linear  in  £j  and  E2).  It  is  known  that  stiffness  decreases  with  fatigue, 
but  detailed  knowledge  of  the  mechanisms  that  explain  this  decrease  is  not  available. 
Therefore,  to  account  for  the  stiffness  reduction,  all  properties  are  reduced  equally. 
This  approximation  is  justified  by  assuming  micromechanical  damage  occurs  in  the 
matrix,  the  matrix  to  fiber  interface,  and  possibly  even  the  fiber  itself.  This  is  also 
justified  because  the  code  does  not  accurately  account  for  delamination  effects,  which 
can  be  a  fatigue  driver  in  some  ply  layups  (35). 

It  is  not  feasible  to  run  the  code  for  100,000  cycles,  as  the  analysis  is  time 
consuming  even  for  a  few  cycles.  Instead,  the  stiffness  reduction  factor  is  scaled  to 
represent  a  large  number  of  actual  fatigue  cycles.  As  seen  in  Figure  2.5,  the  secant 
modulus  (stiffness)  reduces  from  approximately  98.5  percent  to  97.5  percent  in  Stage 
II.  This  is  the  region  of  stable  fatigue  degradation  and  the  largest  portion  of  life  for 
the  composite.  In  this  case,  Stage  II  occurs  over  30,000  cycles.  To  account  for  this 
reduction  in  just  a  few  computational  cycles,  the  material  properties,  which  dictate 
the  overall  stiffness,  are  reduced  by  a  factor,  R.  Representing  the  initial  stiffness  as 
Ki,  the  final  stiffness  as  Kf,  and  the  number  of  computational  cycles  as  n,  R  can  be 
determined  by  solving  the  following  equation: 

KtRn  =  Kf  (3.1) 
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For  Stage  II: 


Ki  =  0.985 K 
Kf  =  0.975 K 

K  =  the  normalizing  (initial)  stiffness 

Equation  (3.1)  is  now  written  as: 

0.985 KRn 
Rn 
R 

For  three  computational  cycles  (n  =  3): 

R  =  0.9966 

or  R  =  99.66%.  In  this  way,  one  computational  cycle  represents  10,000  physical 
cycles.  It  is  now  possible  to  analyze  the  state  of  stress  and  the  state  of  damage  of  a 
structure  after  any  number  of  fatigue  cycles  by  changing  the  variables  n  and  R. 

Without  Figure  2.5,  a  comparison  of  the  initial  stiffness  and  the  secant  stiffness 
at  failure  in  the  equilibrium  curve  can  be  used.  This  is  an  approximation,  however, 
because  of  the  development  of  each  stiffness  matrix.  Remember  the  stiffness  matrix  is 
made  up  of  the  individual  matrices  K,  N\,  and  N2,  where  N\  —  f(q )  and  IV2  =  f(q2)- 
For  the  initial  stiffness,  the  displacements  q  are  zero  and  the  Ni  and  N2  matrices 
are  zero.  Therefore,  there  will  be  difference  between  the  initial  and  final  stiffness 
values.  Both  of  the  above  methods  are  approximations  and  may  not  provide  the 
exact  results  in  the  number  of  cycles  predicted.  Ultimately,  the  stiffness  reduction 
parameter,  R  should  be  varied  to  whatever  value  provides  the  desired  results. 


=  0.975 K 
0.975 
0.985 

=  70.98985 
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The  key  to  accounting  for  damage  progression  and  fatigue  are  the  two  stiff¬ 
ness  reduction  mechanisms:  1)  damage  from  failure  criterion  and  2)  cyclic  loading. 
Damage  from  the  failure  criteria  is  analogous  to  macromechanical  damage,  while 
it  is  assumed  cyclic  loading  produces  micromechanical  damage  leading  to  stiffness 
reduction.  In  actuality,  macromechanical  damage  should  eventually  occur  from  the 
individual  micro-damage. 

3.5  The  Finite  Element  Model 

The  model  chosen  for  this  analysis  is  a  cylindrical  panel  shown  in  Figure  3.5. 
The  sides  measure  30.48  cm  by  30.48  cm  [12  in.  x  12  in.]  and  the  radius  is  also 
30.48  cm  [12  in.].  A  10.16  cm  x  10.16  cm  [4  in.  x  4  in.]  cutout  is  placed  in  the 
middle  of  the  panel.  This  is  the  same  configuration  used  by  Hatfield  (6),  and  allows 
comparison  of  compression  values  with  his  research.  Both  cross-ply  ([0/90];s,  where 
i  =  2  and  6)  and  quasi-isotropic  ([0/45/ —  45/90]jS,  where  i  —  1  and  3)  configurations 
are  used,  resulting  in  layups  of  8  and  24  plies. 

3.5.1  Material  Properties  and  Boundary  Conditions. 

Graphite/epoxy  AS4/3501-6  was  used  for  the  analysis.  The  material  properties  for 
this  composite  are  shown  in  Table  3.1  (6),  (20): 

The  boundary  conditions  were  clamped  along  the  top  and  bottom  edges,  as  this 
resembles  the  experimental  work  of  Hatfield.  This  also  models  the  grips  commonly 
used  in  experimentation,  if  further  testing  is  performed  for  comparison.  Even  though 
the  ends  were  clamped,  movement  in  the  x-direction  only  is  allowed  at  one  end  to 
apply  an  axial  load  or  displacement. 

For  the  static  analysis,  displacement  control  was  used  to  determine  the  maxi¬ 
mum  failure  value.  The  maximum  load  was  determined,  and  fatigue  loads  ranging 
from  30-60%  of  this  value  was  applied  using  load  control.  The  main  reason  for  using 
load  control  is  to  apply  a  constant  stress  to  the  shell,  which  is  the  conventional  S-N 
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Table  3.1  AS4/3501-6  Graphite  Epoxy  Material  Properties 


Property 

SI 

English 

Ei 

135.8  GPa 

19.70xl06  psi 

_  CO 

cq 

II 

cq 

10.89  GPa 

1.579xl06  psi 

Gl2  —  G 13 

6.378  GPa 

9.25xl05  psi 

G23 

3.185  GPa 

4.62xl05  psi 

V\1 

0.276 

0.276 

tply 

0.131  mm 

0.00514  in. 

GpN 

1.5  GPa 

2.176xl05  psi 

&FS 

0.22  GPa 

31.91xl03  psi 

Cm  NT 

43.8  MPa 

6.352xl03  psi 

O'MNC 

43.8  MPa 

6.352xl03  psi 

C  MS 

43.8  MPa 

6.352xl03  psi 

CDN 

50.0  MPa 

7.251xl03  psi 

CDS 

86.0  MPa 

12.473xl03  psi 

€n 

0.011 

0.011 

tr 

0.0065 

0.0065 

approach.  If  displacement  control  were  used,  the  applied  stress  would  decrease  as 
the  structure  became  more  flexible  due  to  damage. 

For  the  given  configuration,  a  uniformly  applied  displacement  will  not  result 
in  a  uniformly  applied  load.  This  is  due  to  the  large  cutout,  where  the  structure 
is  more  flexible  in  this  region.  Therefore,  the  distributed  loads  along  the  edge  for 
displacement  control  will  be  the  greatest  at  the  ends  and  the  least  in  the  middle. 
Likewise  for  load  control,  a  uniform  stress  distribution  along  the  top  edge  will  result 
in  greater  displacements  in  the  middle  of  the  panel,  where  the  cutout  reduces  the 
stiffness.  An  example  is  shown  in  Figure  3.6  of  the  resulting  displacement  from  a 
uniform  applied  edge  load.  Also  plotted  is  the  average  displacement  value.  This 
displacement  difference  can  cause  convergence  difficulty  if  this  distortion  becomes 
too  large. 

It  is  beneficial  to  note  how  a  uniform  stress  is  distributed  to  the  appropriate 
nodes  in  the  FEM  code  for  load  control.  Two  possible  methods  are  consistent  nodal 
loads  and  equivalent  nodal  loads  (41).  Equivalent  nodal  loads  is  a  crude  method 
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Figure  3.6  Actual  and  Average  Displacement  Along  an  Edge  for  a  Uniformly  Ap¬ 
plied  Stress 

where  the  stress  is  converted  into  a  force.  This  force  value  is  then  divided  by  the 
number  of  nodes  within  the  application  region,  resulting  in  the  load  being  divided 
equally  between  all  the  nodes.  There  is  little  or  no  guarantee  that  this  method  will 
produce  a  uniform  stress  field. 

For  the  consistent  nodal  load  method,  the  degrees  of  freedom  are  multiplied 
by  the  shape  functions.  This  product  is  then  integrated  along  the  length  (for  an 
edge  load)  to  produce  the  resulting  consistent  load  value  at  each  node  (41).  For  this 
research,  an  axially  applied  load  will  only  vary  with  the  degree  of  freedom  u.  The 
finite  element  is  eight-noded,  and  the  interpolation  of  u  is  quadratic.  Integrating 
the  quadratic  expression  for  u  along  an  element  edge  will  result  in  a  parabolic  rep¬ 
resentation,  with  the  greatest  load  being  at  the  mid-side  node  (41).  This  parabolic 
representation  has  already  been  programmed  into  the  code.  Since  the  vertex  nodes 
(except  for  the  nodes  at  each  end  of  the  edge)  are  shared  by  two  elements,  their 
individual  load  values  are  summed.  This  results  in  a  consistent  nodal  load  represen- 
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tation  for  the  edge  load  encountered  in  this  research.  This  is  also  true  for  most  of 
the  loading  conditions  for  which  the  FEM  code  was  designed. 

3.5.2  Convergence.  A  desirable  feature  of  finite  element  analysis 
common  to  compatible  elements,  including  this  element,  is  that  as  the  number  of 
elements  increase,  the  solution  converges  to  the  exact  answer.  A  model  with  too  few 
elements  will  present  answers  more  stiff  than  the  actual  structure.  This  is  because 
the  model  is  only  allowed  to  behave  according  to  the  interpolation  functions  chosen 
(41).  To  gain  precision,  either  a  higher-order  interpolation  function  is  required  or 
more  elements  are  needed.  However,  an  increase  in  the  number  of  elements  also 
increases  the  complexity  and  computational  time. 

For  the  compression  case,  the  shell  fails  due  to  geometric  instability.  Hatfield 
(6)  used  Dennis’  (46)  findings  that  an  element  of  the  size  1.27  cm  x  1.27  cm  [0.5 
in.  x  0.5  in.]  provides  sufficient  accuracy  for  this  type  of  problem,  which  is  the  size 
element  also  used  for  the  present  compression  cases.  This  results  in  a  24  element  by 
24  element  array  for  a  panel  with  no  cutout.  The  interior  8x8  elements  are  removed 
for  the  cutout. 

To  find  a  balance  between  convergence  and  complexity  for  the  tension  case, 
a  convergence  analysis  was  performed.  An  8-ply,  cross-ply  configuration  was  used 
for  the  test.  The  convergence  test  started  with  a  6  x  6  array.  As  mentioned  above, 
displacement  control  is  chosen  for  the  static  analysis  to  determine  the  maximum  point 
on  the  equilibrium  curve.  The  number  of  elements  along  the  edges  were  doubled  until 
a  suitable  value  was  found.  The  first  three  arrays  are  plotted  in  Figure  3.7.  Symmetry 
was  used  to  increase  the  number  of  elements.  Symmetry  in  this  manner  only  applies 
to  the  tension  loading,  since  the  compression  loading  fails  by  geometric  instability 
involving  large  out  of  plane  displacements.  For  the  current  boundary  conditions, 
symmetry  exists  about  the  vertical  axis  as  shown  in  Figure  3.8.  Continuity  must 
exist  across  the  arc  for  symmetry  along  this  vertical  axis.  Therefore,  along  the 
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Figure  3.7  Normalized  Equilibrium  Curves  of  Different  Size  Arrays 

vertical  symmetry  line,  the  displacement  v  must  equal  zero.  The  slope  must  also 
be  continuous  across  this  boundary,  so  w,y  =  0  is  required.  By  the  same  token,  the 
rotation  must  also  equal  zero. 

For  symmetry  on  the  horizontal  axis,  at  a;  =  |r,  the  loading  must  also  be 
symmetric.  Instead  of  clamping  one  end  and  loading  the  other,  it  is  necessary  to 
load  both  ends  by  half  the  original  amount.  This  results  in  u  being  zero  at  the 
horizontal  symmetry  axis,  which  is  required.  Continuity  also  requires  that  w,x ,  and 
be  equal  to  zero  on  this  line. 

Before  implementing  quarter-panel  symmetric  models,  it  is  necessary  to  vali¬ 
date  the  above  conditions.  Analysis  of  a  full  24  element  by  24  element  panel  deter¬ 
mines  that  symmetry  does  exist.  The  necessary  null  degrees  of  freedom  along  the 
respective  lines  produce  values  on  the  order  of  10-16,  which  is  assumed  a  negligi¬ 
ble  amount.  The  possibility  exists  that  once  failure  occurs,  this  symmetry  is  lost. 
However,  since  the  panel  exhibits  symmetry,  the  loading  and  hence  the  failure  ap- 
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Figure  3.8  Shell  Symmetry  Axes  under  a  Vertical  Tension  Load 
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Figure  3.9  Comparison  of  a  Full  Panel  and  a  Symmetric  Quarter  Panel 

proximates  symmetry  or  very  close  to  it.  Checking  the  displacement  values  near 
failure  show  that  symmetry  does  begin  to  break  down.  Terms  that  were  formerly  on 
the  order  of  10~16  are  now  on  the  order  of  10-4.  These  terms  are  still  one  to  two 
orders  of  magnitude  smaller  than  the  other  displacements,  though,  and  symmetry  is 
still  assumed  to  exist. 

Symmetry  may  be  verified  by  comparing  a  24  element  by  24  element  cross- 
ply  full  panel  with  a  symmetrically  loaded  quarter  panel  containing  the  same  size 
elements.  The  equilibrium  curves  of  both  panels  is  shown  in  Figure  3.9.  In  order  to 
perform  the  quarter  panel  analysis,  the  boundary  conditions  had  to  be  changed.  The 
null  degrees  of  freedom  were  enforced  along  the  respective  symmetry  axes,  and  the 
appropriate  symmetric  loading  condition  was  applied  to  the  end.  From  the  figure,  it 
is  evident  that  both  panels  behave  almost  identically.  Therefore,  it  is  assumed  that 
a  full  panel  may  be  represented  through  symmetry  by  a  quarter  panel. 
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Figure  3.10  Normalized  Equilibrium  Curves  of  More  Refined  Arrays 

It  is  now  possible  to  continue  the  convergence  check.  The  number  of  elements 
currently  allowed  by  the  program  (without  extensive  modification)  is  600  elements, 
which  will  allow  a  24  by  24  element  array.  A  24  by  24  element  symmetric  quarter 
panel  models  a  48  by  48  element  full  panel  array  with  an  element  size  of  6.35  mm  by 
6.35  mm  [0.25  in.  x  0.25  in.].  Two  quarter  panel  cases  were  chosen  and  are  plotted 
in  Figure  3.10. 

Table  3.2  presents  some  convergence  comparisons.  In  general,  the  table  and 
corresponding  figures  identify  the  trend  that  an  increase  in  the  number  of  elements 
converges  toward  a  less  stiff,  more  exact  answer.  The  ultimate  failure  value  decreases 
with  an  increase  in  the  number  of  elements.  There  are  some  anomalies  in  the  chart, 
however,  such  as  the  small  reduction  in  load  between  the  6  by  6  and  the  12  by  12 
element  arrays.  Another  is  the  large  reduction  is  the  displacement  value  at  which 
failure  occurs.  Therefore,  it  is  advisable  to  re-evaluate  the  convergence. 
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Array  Size 

Normalized 
Load  (%) 

%  Reduction 

Normalized 
Displacement  (%) 

%  Reduction 

6x6 

100.0 

- 

100.0 

- 

12  x  12 

97.6 

2.4 

94.4 

5.6 

24  x  24 

88.5 

9.4 

83.3 

11.8 

30  x  30  (j  panel) 

81.0 

6.3 

83.3 

0.0 

48  x  48  (|  panel) 

77.3 

4.6 

72.2 

13.3 

Table  3.2  A  Comparison  of  Array  Size  with  Normalized  Load  and  Displacement 


Application  of  the  failure  criterion  plays  a  crucial  role  in  determining  the  ul¬ 
timate  panel  failure  load.  In  order  for  the  failure  criterion  to  accurately  represent 
the  actual  material  response,  there  must  be  an  accurate  representation  of  the  stress. 
The  cutout  is  a  large  stress-riser.  From  the  elasticity  approach,  the  stress  approaches 
infinity  near  the  edge  of  the  cutout.  The  elements  should  be  small  enough  to  capture 
an  increase  in  stress  near  the  cutout.  To  analyze  this,  a  local  convergence  test  is 
performed.  The  stress  along  a  ray  from  the  corner  of  the  panel  to  the  corner  of 
the  cutout  is  plotted.  This  is  performed  for  the  various  array  sizes  and  is  shown  in 
Figure  3.11. 

A  singularity  does  exist  at  the  sharp  cutout  corner,  however,  and  the  assump¬ 
tions  made  for  the  finite  element  analysis  begin  to  break  down.  Again,  FEM  must 
meet  the  necessary  requirements  on  an  average  sense,  and  not  necessarily  at  each 
element.  Therefore,  stresses  in  the  element  adjacent  to  the  cutout  will  be  in  error.  In 
fact,  a  large  element  adjacent  to  the  cutout  will  “smear  out”  some  of  the  effect  of  the 
singularity  and  its  resulting  inequality.  A  balance  between  element  size  is  preferred. 

From  Figure  3.11,  it  is  observed  that  the  6x6  and  12  x  12  arrays  actually 
decrease  in  stress  near  the  cutout.  Clearly,  elements  of  this  size  are  too  crude.  The 
48  x  48  element  array  displays  a  large  increase  in  stress  near  the  cutout  corner.  At 
a  distance  of  about  1.5  cm  [0.6  in.]  away  from  the  cutout  corner,  however,  the  48 
x  48,  30  x  30,  and  24  x  24  element  arrays  represent  approximately  the  same  stress 
value. 
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From  this  comparison,  it  is  observed  the  failure  criteria  will  be  more  effective  if 
smaller  elements  are  used.  As  the  element  size  decreases,  a  higher  change  in  stress  can 
be  detected — resulting  in  failure  at  a  lower  global  input  value.  It  is  this  failure  that 
is  driving  down  the  displacement  at  which  ultimate  failure  occurs  from  the  previous 
figures.  For  most  composite  configurations,  the  fiber  is  the  “critical  element”  (28). 
The  residual  strength  of  the  composite  relies  on  this  critical  element.  Having  fibers 
fail  at  a  lower  load — even  in  a  smaller  element — causes  further  element  and  eventually 
ultimate  panel  failure  at  a  lower  load  and  displacement  value.  Therefore,  a  model 
with  smaller  elements  does  a  better  job  of  approximating  an  actual  structure. 

The  argument  against  smaller  elements  is  computer  run-time.  Even  using 
“banded”  matrices,  computer  runs  take  a  long  time  due  to  the  nonlinear  nature  of 
the  analysis  and  the  amount  of  information  that  must  be  generated.  For  a  6  x  6 
element  array,  one  increment  often  takes  10  minutes  on  a  Spark  20  Sun  Workstation. 
If  failure  occurs,  this  increment  must  be  re-applied.  For  a  24  x  24  array,  a  problem 
using  failure  criteria  and  involving  10  increments  often  takes  120  hours  or  more 
to  complete.  In  developing  the  technique  of  handling  fatigue  failure,  it  is  deemed 
necessary  to  reduce  the  size  of  the  model  if  possible.  A  12  x  12  array  has  a  run-time 
of  approximately  24-36  hours.  Using  symmetry  and  a  1.27  cm  x  1.27  cm  [0.5  in.  x 
0.5  in.]  size  element  is  a  good  balance  between  accuracy  and  efficiency.  It  should  be 
sufficient  to  identify  trends  in  the  analysis  in  a  reasonable  amount  of  time. 
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Figure  3.11  Stress  Distribution  for  Various  Array  Sizes 


IV.  Results  and  Discussion 


4.1  Shell  Static  Failure 

Before  the  fatigue  results  are  analyzed,  it  is  important  to  observe  the  results  of 
the  static  loads.  These  will  identify  trends  in  the  panel  behavior  as  well  as  establish 
parameters  for  the  fatigue  loading. 

4.1.1  Static  Compression.  In  compression,  the  shell  fails  from 
instability,  or  collapse.  In  this  failure  mode,  the  shell  deflects  axially  until  it  is  no 
longer  able  to  store  the  applied  energy.  Then,  a  large  out-of-plane  deflection  results 
from  the  axially  applied  load.  It  is  determined  that  very  little  material  failure  occurs, 
although  there  is  some  localized  matrix  failure  due  to  compression.  This  material 
failure  does  little  to  change  the  panel’s  ultimate  collapse  load.  Comparison  with 
Hatfield’s  (6)  computational  results  in  Table  4.1  shows  a  7.5%  reduction  in  collapse 
load,  but  is  still  14  %  above  the  experimental  results  for  the  [0/90]2s  configuration. 
Also,  the  displacement  at  collapse  has  a  slight  but  negative  change  from  computa¬ 
tional  results  that  do  not  consider  material  failure.  The  quasi-isotropic  ply  yields 
similar  results. 


Table  4.1  The  Effect  of  Material  Failure  in  Collapse  Analysis 


Type 

Collapse  Load 
(N)  [lbs] 

Displacement  at  Collapse 
(mm)  [in.] 

Cross-ply  with  Material  Failure 

5,831  [1,311] 

0.25  [0.010] 

Cross-ply  without  Material  Failure  (6) 

6,307  [1,418] 

0.28  [0.011] 

Cross-ply  Experimental  Results  (6) 

5,115  [1,150] 

1.02  [0.040] 

Quasi-isotropic  with  Material  Failure 

5,747  [1,292] 

0.22  [0.009] 

Quasi-isotropic  without  Material  Failure  (6) 

6,272  [1,410] 

0.25  [0.010] 

Quasi-isotropic  Experimental  Results  (6) 

4,403  [990] 

0.711  [0.028] 

4.1.2  Static  Tension.  Material  failure  is  the  primary  cause  of 
panel  failure  for  the  tension  case.  This  is  displayed  through  the  determination  of 
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Figure  4.1  FEM  Model  of  Plate  used  for  Comparison  with  Jen,  et  al. 


Configuration 

Jen,  et  al. 
(MPa)  [ksi] 

FEM  Analysis 
(MPa)  [ksi] 

%  Difference 

[90/ 0]4s 

258  [37.4] 

260  [37.7] 

0.8 

[±45]4s 

73  [10.6] 

106  [15.4] 

45.4 

[0/90/  ±  45]2s 

203  [29.4] 

225  [32.6] 

10.8 

[90/0/  ±  45]2s 

141  [20.4] 

225  [32.6] 

59.6 

Table  4.2  Comparison  of  Results  for  Flat  Plate  with  a  Hole 


ultimate  load,  displacement  at  ultimate  load,  and  failure  progression  through  the 
panel  for  the  static  case. 

In  order  to  validate  the  approach,  comparison  is  made  with  experimental  re¬ 
sults.  Results  for  the  curved  shell  in  tension  were  not  available,  so  comparison  is 
made  with  experimental  results  for  a  flat  plate  from  Jen,  et  al.  (34).  A  12.0  cm  x 
2.75  cm  [4.72  in.  x  1.08  in.]  plate  with  a  0.4  cm  [0.16  in.]  diameter  hole  was  com¬ 
pared  for  4  different  configurations.  Symmetry  was  used  to  develop  a  quarter-plate 
FEM  model,  which  is  shown  in  Figure  4.1.  The  comparisons  are  in  Table  4.2. 

The  [90 / 0]4S  configuration  is  extremely  close  to  the  actual  value,  and  this  de¬ 
velops  confidence  in  the  ability  of  the  Hashin  failure  criteria  to  capture  the  damage 
within  the  plate.  There  is  a  much  larger  variance  in  the  results  for  the  [+45/  —  45]4s 
and  [90/0/ ±45]^  configurations.  This  is  because  delamination  plays  a  larger  role  in 
the  failure  of  these  configurations.  As  mentioned  in  Chapter  3,  the  current  approach 
does  not  capture  delamination  failure  well. 

Just  by  varying  the  ply  layups  in  the  quasi-isotropic  configurations  causes  an 
increase  in  delamination  initiation.  Jen  (34)  states  that  delamination  begins  for  the 
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Table  4.3  Static  Analysis  Results  for  8-ply  Graphite/Epoxy  under  Axial  Tension 


Configuration 

Ultimate  Load 
(N)  [lbs] 

Displacement  at  Ultimate  Load 
(mm)  [in.] 

Cross-ply 

124,958  [28,092] 

2.54  [0.10] 

Quasi-isotropic 

99,645  [22,401] 

3.048  [0.12] 

[0/90/  ±  45]2s  case  between  the  0°  and  90°  plies,  but  almost  at  the  panel  ultimate 
stress.  For  the  [90/0/±45]2s  case,  delamination  begins  between  the  0°  and  +45°  plies. 
This  is  because  the  0°  ply  is  the  primary  axial  load  carrying  ply.  Having  the  +45° 
ply  adjacent  to  this  ply  causes  a  large  shear  stress  between  the  two  plies,  resulting 
in  delamination  initiation  at  a  lower  stress.  For  cases  where  delamination  is  not  the 
primary  failure  mode  of  the  composite,  the  FEM  code  with  the  Hashin  criterion  does 
a  good  job  of  modeling  the  failure  progression  and  capturing  the  ultimate  stress. 

An  8-ply  cross-ply  ([0/90]2s)  and  a  quasi-isotropic  ([0/  ±  45/90]s)  shell  were 
analyzed  quasi-statically  until  failure.  Displacement  control  was  used  to  determine 
the  panel  ultimate  load  value,  and  the  results  are  listed  in  Table  4.3.  The  equilibrium 
curves  for  these  configurations  are  shown  in  Figure  4.2. 

The  progression  of  failure  through  a  cross-ply  panel  is  plotted  in  Figures  4.3 
and  4.4.  For  certain  points  along  the  equilibrium  curves  shown  in  Figure  4.2,  a 
“snapshot”  of  failure  is  taken.  The  data  shows  that  membrane  stresses  dominate 
over  bending  stresses  in  tension  loading  in  the  panel.  That  is  to  say,  for  a  certain 
ply  orientation  the  stress  is  approximately  uniform  through  the  thickness.  Failure 
is  not  distinguished  for  individual  plies,  but  plies  of  one  orientation  tend  to  show 
the  same  amount  of  failure.  Therefore,  if  matrix  failure-  occurs  in  a  90°  ply,  it  often 
occurs  in  every  90°  ply.  The  same  applies  to  fiber  failure  in  the  0°  plies.  There  are  a 
few  exceptions  to  this,  but  they  are  usually  localized  and  away  from  the  progressing 
damage  front.  Also,  failure  is  determined  to  exist  in  the  ply  if  at  least  50%  of  the 
Gauss  points  within  that  ply  have  failed. 
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Figure  4.2  Equilibrium  Curves  for  Two  Laminates  under  Axial  Tension 

Failure  in  the  0°  plies  is  plotted  in  Figure  4.3.  The  primary  failure  mode  is  fiber 
failure,  but  the  matrix  may  also  fail  due  to  tension  in  these  plies.  Figure  (a)  shows 
that  first  fiber  failure  occurs  near  the  cutout  as  expected,  and  at  a  very  high  load. 
Fiber  failure  extends  vertically  along  the  edge  of  the  cutout,  as  well  as  progressing 
across  the  width  in  (b),  which  is  perpendicular  to  the  primary  load  direction.  Also 
notice  some  matrix  failure  due  to  compression  occurs  along  the  “hidden”  edge  of 
the  cutout.  Figure  (b)  is  at  the  panel’s  ultimate  or  maximum  load.  As  shown  in 
figures  (c)  and  (d),  just  a  small  increase  in  displacement  causes  the  damage  front  to 
progress  across  the  panel,  leading  to  panel  failure. 

For  the  90°  plies,  failure  is  already  evident  at  a  displacement  of  0.635  mm  [0.025 
in.],  corresponding  to  an  axially  applied  load  of  34.2  kN  [7,683  lbs.].  As  shown  in 
Figure  4.4(a),  this  is  a  matrix  failure  due  to  tension.  Once  the  displacement  reaches 
1.143  mm  [0.045  in.],  the  characteristic  damage  state  (CDS)  has  occurred,  as  shown 
in  (b).  Matrix  cracking  has  occurred  and  stabilized  in  the  panel.  A  further  increase 
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in  applied  displacement  only  results  in  a  few  more  matrix  failures  for  the  90°  plies, 
as  seen  in  figure  (c). 

As  mentioned  in  Chapter  1,  there  are  different  approaches  to  analyzing  failure 
in  a  composite.  One  approach  is  the  fracture  mechanics  approach,  which  bases 
failure  on  a  crack  type  of  damage.  While  Figure  4.3  does  shown  a  crack-like  damage 
progression,  the  damage  occurs  over  the  width  of  several  elements.  This  leads  to 
another  approach  for  failure,  the  damage  growth  method.  Because  of  the  different 
constituents  and  orientations  in  composite  laminates,  failure  is  often  spread  across 
an  area.  This  is  especially  obvious  for  the  90°  ply  failures  in  Figure  4.4.  While  this 
is  not  the  primary  damage  leading  to  panel  failure,  the  combination  of  the  fiber  and 
matrix  failures  over  several  elements  validates  use  of  a  damage  area  or  accumulation 
approach. 

Similar  results  are  plotted  for  a  quasi-isotropic  8- ply  shell  in  Figures  4.5  -  4.8. 
From  these  figures,  it  is  seen  that  the  damage  pattern  takes  a  different  shape.  This 
configuration  lends  itself  well  to  the  damage  progression  failure  model.  Again,  failure 
originates  in  the  90°  plies  near  the  cutout.  The  —45°  plies  are  the  next  to  show  large- 
scale  failure.  These  plies  are  adjacent  to  the  90°  plies,  and  must  pick  up  the  stress 
distribution  due  to  the  early  90°  ply  matrix  failure.  The  +45°  plies  fail  in  matrix 
next,  and  finally,  the  fibers  in  the  0°  plies  fail.  The  failure  pattern  spreads  out  in  a 
45°  pattern,  including  the  fiber  failures.  In  general,  this  damage  pattern  correlates 
with  experimental  failure  patterns,  such  as  that  found  in  Daniels,  et  al.  (47). 

It  is  important  to  note  that  in  both  configurations,  delamination  was  not  shown 
to  occur.  According  to  Barboni  et  al.(35),  delamination  is  not  very  prevalent  in  cross- 
ply  laminates.  However,  delamination  did  not  occur  for  the  quasi-isotropic  panel 
either,  nor  after  a  large  amount  of  failure  occurred  in  either  the  panel.  This  departs 
from  the  common  knowledge  that  delamination  is  an  active  damage  mechanism  in 
laminate  composites. 
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(a)  Displ.=l. 905mm,  Load=95.5kN  (b)  Displ.=2.54mm,  Load=125.0kN 


(c)  Displ.=2.794mm,  Load=115.3kN  (d)  Displ.=3.048mm,  Load=119.7kN 

Figure  4.3  Failure  Progression  for  0°  Plies  in  a  Cross-ply  Shell 
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(a)  Displ.— 0.635mm,  Load=34.2kN 


(b)  Displ.=l. 143mm,  Load=61.0kN 


(c)  Displ. =1. 905mm,  Load=95.5kN 

Figure  4.4  Failure  Progression  for  90°  Plies  in  a  Cross-ply  Shell 


HI 


(a)  Displ.=l. 905mm,  Load=74.7kN  (b)  Displ.=2.54mm,  Load=91.4kN 

Figure  4.7  Failure  Progression  for  —45°  Plies  in  a  Quasi-Isotropic  Shell 

There  are  several  reasons  why  delamination  did  not  appear  in  the  results,  which 
have  been  mentioned  previously.  From  the  chosen  delamination  failure  criterion,  the 
stresses  contributing  to  delamination  are  <733,023,  and  <7i3-  Again,  033  is  assumed  to 
be  zero  by  the  theory.  Therefore,  only  a23  and  013  (04  and  a5  in  our  notation)  are 
capable  of  causing  delamination.  If  free-edge  effects  were  considered,  these  stresses 
would  increase  near  the  boundaries.  Free-edge  effects  are  not  considered  in  the 
theory,  so  this  increased  stress  distribution  due  to  the  free-edge  is  not  captured. 
Lee  (18)  also  encountered  this  problem  with  the  criterion  he  developed,  stating  an 
inability  to  refine  the  mesh  near  the  edge  to  capture  these  effects. 


4-2  Shell  Fatigue  Loading 

Load  control  is  used  in  order  to  apply  a  constant  stress.  After  the  first  cycle, 
the  stiffness  is  arbitrarily  reduced  to  model  experimental  fatigue  results.  Load  con¬ 
trol  results  in  an  uneven  displacement  along  the  edge  where  it  is  applied.  This  is  due 
to  the  reduced  panel  stiffness  in  the  area  of  the  cutout.  In  order  to  plot  the  equilib- 
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(a)  Displ.=0. 635mm,  Load=25.9kN 


(b)  Displ.=1.27mm,  Load=51.4kN 


Matrix- T  ension 
HI  Fiber-Tension 

Matrix-Compression 


(c)  Displ.=l. 905mm,  Load=74.7kN 

Figure  4.8  Failure  Progression  for  90°  Plies  in  a  Quasi-Isotropic  Shell 

rium  curve,  the  average  displacement  across  the  edge  is  used.  Several  combinations 
of  load  value,  stiffness  reduction  parameter,  number  of  cycles,  and  ply  thicknesses 
were  evaluated.  This  is  a  benefit  of  this  approach,  since  the  parameters  may  be 
adjusted  to  highlight  certain  effects.  Due  to  similarity,  only  the  results  for  the  8-ply 
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Table  4.4  Fatigue  Parameters  for  the  Compression  Equilibrium  Curves 


Configuration 

Fatigue  Load  as  %  of 
Collapse  Load 

Stiffness  Reduction 
Parameter,  R 

Number  of  Computer 
Cycles  Applied 

Cross-ply 

29 

0.90 

6 

Quasi-isotropic 

29 

0.995 

6 

configurations  are  shown — for  both  the  cross-ply  and  the  quasi-isotropic  cases.  Equi¬ 
librium  curves  for  compression  and  tension  are  plotted  for  a  representative  model 
for  each  configuration. 

4-2.1  Compression.  In  compression,  panel  failure  is  mainly  due 
to  geometric  nonlinearities  and  not  material  failure.  The  compression  loadings  in 
fatigue  were  very  sensitive  to  the  fatigue  stiffness  reduction  parameter,  R.  If  R 
were  too  large,  the  FEM  code  would  not  converge,  even  on  the  first  fatigue  cycle. 
If  it  were  too  small,  final  panel  failure  due  to  fatigue  would  not  occur.  Table  4.4 
lists  the  fatigue  parameters  for  the  configurations  corresponding  with  the  displayed 
equilibrium  curves. 

The  effects  of  the  stiffness  reduction  can  be  seen  in  the  increase  in  axial  dis¬ 
placement  that  occurs  each  computer  cycle  in  Figures  4.9  and  4.10.  It  is  important 
to  remember  that  one  computer  cycle  represents  thousands  of  actual  fatigue  load¬ 
ing  cycles.  For  the  cross- ply  panel,  the  larger  change  in  stiffness  each  cycle  shows 
a  large  change  in  axial  displacement.  The  stiffness  reduction  for  this  chart  is  0.90. 
The  maximum  axial  displacement  reached  is  0.082  mm  [0.003  in.].  The  displacement 
at  collapse  for  this  panel  is  0.25  mm  [0.009  in.],  so  many  more  cycles  are  required 
before  collapse  occurs  by  fatigue  for  this  applied  load.  Comparing  the  change  in 
stiffness  values  each  cycle  shows  a  90%  reduction,  which  correlates  exactly  with  the 
stiffness  reduction  parameter.  Since  the  reduction  parameter,  R,  is  0.995  for  the 
quasi-isotropic  panel,  there  is  very  little  change  in  displacement  each  cycle. 

After  these  runs,  the  conditions  were  changed  to  show  final  collapse  in  fatigue. 
This  was  achieved  in  four  cycles  for  both  configuration  at  a  fatigue  load  60%  of  the 
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Figure  4.9  Fatigue-Compression  Cycles  for  a  Cross-ply  Shell- (29%  Load) 


Figure  4.10  Fatigue-Compression  Cycles  for  a  Quasi-isotropic  Shell-(29%  Load) 
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Figure  4.11  Fatigue-Compression  Cycles  for  a  Cross-ply  Shell- (60%  Load) 

collapse  load  and  an  R  of  0.88.  For  the  cross-ply  configuration,  the  displacement  at 
collapse  is  0.15  mm,  shown  in  Figure  4.11.  This  is  much  lower  than  the  quasi-static 
results  for  displacement  at  collapse,  which  was  0.25  mm.  The  same  effect  occurred 
for  the  quasi-isotropic  panel,  so  those  results  are  not  shown.  The  fatigue  collapse 
displacement  was  at  0.17  mm,  compared  to  the  static  collapse  displacement  of  0.22 
mm. 


4.2.2  Tension.  In  tension,  the  individual  failures  cause  a  change  in 
slope  in  the  equilibrium  curve.  This  material  nonlinearity  can  be  accounted  for  by 
applying  the  failure  criterion  and  iterating  until  the  new  equilibrium  value  is  reached. 
To  account  for  the  stiffness  reduction  due  to  fatigue,  the  stiffness  terms  are  reduced 
by  a  set  amount  each  cycle.  Again,  several  configurations  and  parameter  values  are 
used,  and  the  results  for  the  8-ply  cross- ply  and  quasi-isotropic  panels  are  displayed 
to  show  the  trends.  The  fatigue  parameters  for  the  displayed  configurations  are  listed 
in  Table  4.5. 
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Table  4.5 

Fatigue  Parameters  for  the  Tension  Equilil 

brium  Curves 

Configuration 

Fatigue  Load  as  %  of 
Ultimate  Load 

Stiffness  Reduction 
Parameter,  R 

Number  of  Computer 
Cycles  Applied 

Cross-ply 

30 

0.88 

11 

Quasi-isotropic 

32 

0.88 

11 

Figures  4.12  and  4.13  display  the  different  equilibrium  curves  for  the  tension 
case.  Again,  cycle  is  referring  to  a  computer  load  cycle  which  represents  thousands 
of  actual  fatigue  cycles.  From  these  curves,  it  is  observed  that  panel  failure  is  never 
reached.  Even  though  the  applied  load  is  low  (30%  of  ultimate),  it  is  high  enough 
to  generate  localized  failure  near  the  cutout  for  the  first  cycle.  If  a  higher  load  were 
applied,  or  a  larger  stiffness  reduction  introduced,  panel  failure  would  occur.  For 
example,  the  stiffness  reduction  parameter,  R,  for  the  cross-ply  panel  was  reduced 
to  0.60,  and  panel  failure  occurred  in  four  computer  cycles.  Since  load  control  is 
employed,  panel  failure  results  in  either  the  convergence  criteria  being  exceeded  or 
a  sudden  occurrence  of  a  very  large  displacement.  This  is  because  the  load  control 
method  cannot  traverse  the  maximum  point  on  the  equilibrium  curve. 

Another  interesting  effect  of  the  fatigue  stiffness  reduction  technique  for  tension 
loading  is  that  no  further  material  failure  occurs.  From  before,  F  =  kd  or  R  =  Kq. 
Since  R  is  applied  at  the  same  value  and  K  is  reduced,  q  increases.  The  stress  in 
the  load  carrying  elements  remains  approximately  the  same.  Therefore,  this  tech¬ 
nique  will  never  generate  macromechanical  critical  failure  for  a  stress-based  failure 
criterion.  An  interesting  result  does  occur,  however.  In  the  “hidden”  region  below 
the  cutout,  no  global  tension  load  is  carried.  The  load  at  this  point  has  “sheared” 
out  around  the  cutout.  In  fact,  results  from  the  static  test  show  compression  in  this 
area.  For  the  analysis,  the  only  material  failure  with  the  Hashin  criterion  for  the 
panel  due  to  fatigue  is  in  this  area.  The  panel’s  global  displacement  is  causing  an 
increased  compressive  stress  in  this  local  region,  which  results  in  matrix  compressive 
failure.  This  region  is  not  the  primary  load-carrying  region,  however,  and  it  does 
not  affect  the  panel’s  ultimate  response. 
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Figure  4.12  Equilibrium  Curves  for  a  Cross-ply  under  Tension-Fatigue 


Figure  4.13  Equilibrium  Curves  for  a  Quasi- isotropic  Panel  under  Tension-Fatigue 
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Since  the  displacement  and  strains  are  increasing  in  each  element,  but  not 
the  stress,  a  strain-based  failure  criterion  was  added.  The  maximum  strain  failure 
criterion  was  used  to  account  for  the  macro-damage  that  leads  to  the  end  of  the 
fatigue  life.  The  first  cycle  is  applied  as  before,  and  the  Hashin  failure  criterion  is 
used.  The  new  equilibrium  point  is  determined  from  this  failure,  and  then  the  load 
is  relaxed.  The  maximum  strain  criterion  is  then  applied  on  the  next  cycle.  To  avoid 
redundancy,  the  strain  criterion  is  not  applied  to  Gauss  points  which  have  already 
failed  (by  the  Hashin  criterion). 

The  8-ply  cross-ply  configuration  was  used  to  test  the  addition  of  the  maximum 
strain  failure  criterion.  The  equilibrium  curves  are  shown  in  Figure  4.14.  For  this 
case,  the  stiffness  reduction  ratio,  R,  was  kept  at  0.88,  and  the  load  was  30%  of  the 
ultimate  panel  load.  With  the  maximum  strain  criterion,  panel  failure  occurred  in 
five  cycles.  Therefore,  this  technique  is  able  to  model  all  three  phases  of  the  fatigue 
life  in  the  axial  tension  case.  As  determined  before,  the  panel  cannot  withstand  much 
fiber  failure.  Fiber  failure  in  strain  occurs  in  the  fifth  cycle.  This  failure  originates  at 
the  element  adjacent  to  the  cutout,  just  as  it  did  in  the  static  case.  Further  element 
failure  in  this  cycle  propagates  out  from  this  element,  resulting  in  element  failure 
from  both  the  Hashin  and  strain  criteria.  In  this  way,  the  maximum  strain  criterion 
is  a  “catalyst”  to  start  the  final  fatigue  failure  present  in  Stage  III. 

Reducing  the  stiffness  of  the  panel  simulates  the  affect  of  fatigue.  The  stress 
distribution  during  the  fatigue  cycle,  after  CDS,  is  determined,  as  well  as  the  panel 
displacement  values.  In  the  fatigue  cycles,  the  displacements  and  strains  are  provided 
at  virtually  any  location.  For  the  compressive  fatigue  loads,  where  the  mode  of  failure 
is  primarily  due  to  geometric  instability,  no  additional  failure  criteria  are  needed.  For 
the  tension  fatigue  loads,  where  material  failure  dominates  the  ultimate  failure  mode, 
the  maximum  strain  criterion  is  used  to  initiate  failure  in  the  critical  constituent- 
which  is  the  0°  fibers  near  the  cutout.  This  technique  is  now  able  to  capture  all  three 
stages  in  the  life  of  a  composite  material  for  either  compression  or  tension. 
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Figure  4.14  Equilibrium  Curves  for  a  Cross-Ply  Panel  in  Fatigue  with  Max  Strain 
Criterion 


V.  Conclusions 

This  research  developed  a  different  approach  to  analyzing  the  stress  of  a  composite 
shell  in  fatigue.  This  approach  was  built  upon  several  proven  concepts,  such  as  the 
SLR  theory  and  the  Hashin  failure  criterion.  Some  key  conclusions  can  be  drawn 
from  this  research: 

1.  Finite  Element  Methods  (FEM)  are  a  useful  tool  to  analyze  composites  due  to 
their  ability  to  determine  the  displacement,  strain,  and  stress  at  virtually  any 
point  in  the  structure  at  any  time  during  loading. 

2.  Failure  criteria  work  well  with  FEM  because  the  stresses  and  strains  are  readily 
available.  A  favorable  criterion  for  composite  laminates  is  the  Hashin  failure 
criterion  because  it  analyzes  multiple  failure  modes  individually.  A  benefit  of 
this  criterion  is  the  ability  to  tailor  stiffness  reduction  to  the  type  of  failure 
that  occurs,  allowing  a  more  accurate  representation  of  the  actual  structure. 
This  does  not  limit  analysis  to  one  type  of  failure  mode,  where  other  failure 
modes  must  be  ignored. 

3.  The  SLR  theory,  along  with  an  algorithm  for  “iterating  on  damage”  allows 
both  material  and  geometric  nonlinearities  to  be  evaluated. 

4.  Even  though  the  code  accounts  for  delamination  failure,  delamination  did  not 
appear  in  any  of  the  failure  methods.  If  a  more  accurate  determination  of 
delamination  is  required,  a  theory  which  accounts  for  free-edge  effects  as  well 
as  stress  normal  to  the  laminate  is  required. 

5.  Fatigue  causes  stiffness  reduction  in  composite  materials.  Stiffness  reduction 
due  to  fatigue  can  be  incorporated  into  existing  FEM  code  to  determine  the 
displacement  and  strain  in  a  composite  cylindrical  shell  with  a  square  cutout. 

6.  By  varying  the  parameters  R  and  n,  it  is  possible  to  match  the  stiffness  re¬ 
duction  of  any  composite  material  that  displays  a  reduction  in  stiffness  due  to 
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fatigue.  This  also  allows  analysis  to  be  tailored  to  any  region  in  the  fatigue 
life  at  any  stress.  Therefore,  the  given  approach  may  be  applied  to  a  myriad 
of  materials,  geometries,  and  load  conditions. 

7.  A  strain-based  criterion  is  able  to  determine  when  panel  failure  occurs  due  to 
a  fatigue  load  in  tension,  by  initiating  macro-damage  in  the  panel. 

8.  This  technique  is  not  limited  to  the  first  ply  failure  or  just  the  initiation  of 
damage,  but  is  also  applicable  after  the  characteristic  damage  state  and  up  to 
the  end  of  the  fatigue  life. 

Further  research  is  necessary  to  validate  this  approach.  It  would  also  be  bene¬ 
ficial  to  perform  experimental  analysis  for  comparison,  as  well  as  having  a  graphical 
post-processor  to  speed  analysis  of  the  results.  However,  this  research  does  present 
a  way  to  analyze  composite  materials  with  complex  geometries  under  fatigue  loads 
and  determines  the  resulting  trends.  The  research  does  not  limit  the  type  of  failure 
to  one  mode,  but  presents  more  of  a  complete  composite  analysis.  Not  only  is  the 
initial  damage  determined,  but  fatigue  effects  are  also  incorporated  as  the  structure 
is  loaded. 
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Appendix  A.  SLR  program  guide, 

version  fshelljb.f 


The  following  is  a  description  of  how  to  compile  and  run  the  program,  the  appropriate 
cards  for  the  input  file,  and  a  sample  input  file.  To  run  the  program,  the  FORTRAN 
file  fshelljb.f  is  compiled  using  a  SUN  FORTRAN  77  compiler.  In  order  for  it  to 
compile,  open  a  terminal  to  the  directory  containing  the  fshelljb.f  file.  It  is  necessary 
for  the  compiler  to  read  an  extended  line  of  up  to  134  characters.  This  is  done 
through  the  extend  command,  -e.  It  is  also  necessary  for  the  compiler  to  read  line 
continuations  well  into  the  300’s  (for  the  stiffness  arrays).  The  command  is  -N1999 
(N,  1  as  in  line,  999  for  up  to  999  continuation  lines).  The  resulting  input  line  is: 

f77  -e  -N1999  fshelljb.f 

Once  the  program  compiles,  the  execution  file  is  “a.out”.  Typing  “a.out”  will 
result  in  a  prompt  line: 

WHAT  IS  THE  INPUT  FILE? 

Type  in  the  input  file  name.  The  input  file  should  be  a  text  file,  but  with  no 
extensions.  This  will  create  two  new  files,  with  the  same  name  but  the  extensions 
“,ans”  and  “.out” .  If  the  nodal  forces  are  to  be  displayed,  they  will  be  displayed  to 
the  screen  as  it  runs  as  well  as  to  the  .ans  file.  The  .out  file  contains  the  input  data 
as  well  as  the  displacements  and  stresses  (if  calculated) . 

The  input  file  is  made  up  of  “cards”  of  data.  Unless  noted,  if  multiple  entries 
are  required  for  a  card,  separate  each  entry  by  a  comma.  Units  are  not  distinguished, 
so  it  is  up  to  the  user  to  ensure  all  input  properties  are  in  compatible  units.  The 
following  key  is  used: 

•  CARD  1:  title 

The  first  line  of  the  program  is  any  text  necessary  to  identify  the  file.  The 
line  may  not  exceed  133  characters. 

•  CARD  2:  iel,  npe,  nanal(l),  nanal(2),  nanal(3),  imesh,  nprnt,  nprint, 
ncut 

iel  is  the  structure  type:  1  for  a  plate  or  beam,  2  for  circular  cylindrical 
shell. 

npe  is  the  nodes  per  element:  either  4  or  8  (NOTE:  Only  an  8-noded 
element  should  be  used  for  a  cylindrical  shell  to  maintain  continuity  between 
elements) . 

nanal(l)  determines  the  type  of  analysis:  0  for  nonlinear  analysis,  1  for 
linear  analysis,  and  2  for  eigenvalue  (bifurcation)  analysis. 
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nanal(2)  is  for  material  type:  0  is  an  arbitrary  laminate,  1  is  for  isotropic 
materials,  2  is  for  symmetric  laminates,  and  3  is  for  sandwich  composite  ma¬ 
terials. 

nanal(3)  is  the  type  of  theory  used:  0  for  SLR  (nonlinear  shell  theory),  1 
for  Von  Karman  plate/Donnell  shell  (linear  theory). 

imesh  mesh  generation:  0  for  manual,  1  for  automatic  (always  the  recom¬ 
mended  selection). 

nprnt  determines  if  the  elasticity  matrices  are  printed:  0  for  no,  1  for  yes. 
Only  the  initial  matrices  are  printed. 

nprint  determines  if  the  element  stiffness  matrices  are  printed:  0  for  no, 
1  for  yes. 

ncut  determines  the  number  of  elements  which  will  be  removed  from  the 
mesh  arrangement:  enter  0  if  no  elements  are  cutout. 

•  CARD  3: 

If  nanal(l)=0  (nonlinear  analysis)  input  intyp,  nine,  imax,  ires,  tol,  fa¬ 
tigue 

intyp  increment  type:  0  for  load  control,  1  for  displacement  control, 
nine  is  the  number  of  increments  of  type  intyp. 

imax  is  the  maximum  number  of  iterations  before  it  is  assumed  conver¬ 
gence  will  not  be  met  (often  on  the  order  of  100) 

ires  determines  how  often  to  update  the  stiffness  matrix:  0  for  every  iter¬ 
ation,  1  for  every  increment 

tol  tolerance  to  be  met  for  convergence,  often  0.1%  or  0.001. 

fatigue  if  fatigue  cycles  with  stiffness  reduction  is  to  be  used,  0  for  no,  1 
for  yes. 

If  nanal(l)=l  (linear  analysis)  skip  this  card. 

If  nanal(l)=2  (eigenanalysis)  read  rstep,  which  is  the  eigenvalue  step. 

•  CARD  4: 

If  nanal(l)=0  (nonlinear  analysis)  and  fatigue=0  list  table  (nine),  which 
is  the  table  of  multipliers  to  determine  each  increment  step  for  either  load  or 
displacement  control. 

If  nanal(l)=0  (nonlinear  analysis)  and  fatigue=l  list  minval,  maxval, 
factor 

minval  is  the  multiplier  for  the  minimum  value  in  the  fatigue  cycle 
maxval  is  the  multiplier  for  the  maximum  value  in  the  fatigue  cycle 
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factor  is  the  factor  that  is  multiplied  to  reduce  the  stiffness  terms  each 
fatigue  cycle.  For  example:  For  a  5%  reduction  in  material  properties  each 
fatigue  cycle,  factor=0.95. 

If  nanal(l)  =  l  (linear  analysis)  skip  this  card. 

If  nanal(l)=2  (eigenanalysis)  skip  this  card. 

•  CARD  5a: 

If  imesh=0  (manual  mesh  generation)  read  in  nem,  nnm,  nx,  ny 
nem  is  the  number  of  elements, 
nnm  is  the  number  of  nodes  in  the  mesh, 
nx  is  the  number  of  elements  in  the  x-direction 
ny  is  the  number  of  elements  in  the  y-  or  s-direction 
If  imesh=l  (automatic  mesh  generation)  just  read  in  nx,  ny 
nx  is  the  number  of  elements  in  the  x-direction 
ny  is  the  number  of  elements  in  the  y-  or  s-direction 

•  CARD  5b: 

If  imesh=0  (manual  mesh  generation)  read  in  nod(i),  x(i),  y(i) 
nod(i)  is  the  node  number 
x(i)  is  nod(i)’s  x  position 
y(i)  is  nod(i)’s  y  or  s  position. 

If  imesh=l  (automatic  mesh  generation)  read  in  dx(i)  and  then  dy(i) 

dx(i)  is  the  array  of  distances  between  the  x  nodes.  There  should  be  nx 
number  of  entries. 

dy(i)  is  the  array  of  distances  between  the  y  nodes.  There  should  be  ny 
number  of  entries. 

•  CARD  6: 

If  ncut=0  (no  cutout  elements)  skip  this  card 

If  ncut=‘n’  (‘n’  number  of  cutout  elements)  read  in  icut(ncut),  which  are 
the  list  of  element  numbers  to  be  cutout.  There  should  be  a  total  of  ncut 
elements  listed. 

•  CARD  7:  LD,PO 

LD  is  the  type  of  load  applied:  0  for  no  load,  1  for  transverse  load,  2  for 
dead  weight,  and  3  for  axial  loading 

PO  is  the  stress  value  applied  for  distributed  loading 
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•  CARD  8:  Unless  LD=3  (axial  loading),  skip  all  of  CARD  8,  including 

CARD  8a  and  8b. 

•  CARD  8a:  For  LD=3  (axial  loading)  read  nedge  which  is  the  number  of 
nodes  with  in-plane  loading 

•  CARD  8b:  For  LD=3  (axial  loading)  read  iedge(nedge),  which  is  the  array 
of  nodes  where  in-plane  loading  occurs. 

•  CARD  9a:  nbdy  is  the  number  of  nodes  where  degrees  of  freedom  are  spec¬ 
ified. 

•  CARD  9b:  nbound  describes  which  degrees  of  freedom  at  which  nodes  are 

specified.  Each  line  of  this  card  contains  eight  entries.  The  first  entry  is  the 
node  number,  and  the  next  seven  entries  correspond  to  the  degrees  of  freedom 
[u,v,w,w,x  ,w,s  If  the  degree  of  freedom  is  not  specified,  a  zero  is 

placed  in  its  space.  If  it  is  to  be  specified,  a  one  is  placed  in  its  corresponding 
position.  For  example,  if  degrees  of  freedom  u,v,  and  w  are  to  be  specified  for 
node  112  and  degrees  of  freedom  w,x  and  w,s  are  to  be  specified  for  node  114, 
the  entries  would  appear: 

112,1,1,1,0,0,0,0 

114,0,0,0,1,1,0,0 

There  should  be  nbdy  number  of  lines  to  this  card.  A  carriage  return  is 
necessary  between  node  entries,  but  a  comma  is  not. 

•  CARD  9c:  vbdy(i)  is  an  array  of  the  prescribed  displacements  from  above. 
They  are  sequentially  placed.  From  the  example  above,  vbdy(i)  would  contain 
five  entries. 

•  CARD  10:  nbsf  is  the  number  of  degrees  of  freedom  containing  arbitrary 
loads. 

•  CARD  10a:  ibsf(i)  is  an  array  of  nbsf  node  numbers  with  prescribed  loads. 
This  card  is  skipped  if  nbsf  is  zero. 

•  CARD  10b:  vbsf(i)  is  an  array  of  the  loads  corresponding  to  the  degrees  of 
freedom  in  ibsf(i),  therefore  sequence  is  important.  Again,  this  card  is  skipped 
if  nbsf  is  zero. 

•  CARD  11:  is  for  the  material  properties. 

If  nanal(2)=l  (isotropic),  enter  ey,  pnu,  ht. 
ey  is  Young’s  modulus, 
pnu  is  Poisson’s  ratio, 
ht  is  the  thickness. 

If  nanal(2)=0,2  (composite  laminate),  enter  el,e2,gl2,pnul2,gl3,g23 
el  is  the  modulus  in  the  1  or  fiber  direction. 
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e2  is  the  2  or  matrix  modulus. 

gl2  is  the  shear  modulus  in  the  1-2  direction. 

pnul2  is  the  Poisson’s  ratio  in  the  1-2  direction 

gl3,  g23  are  the  shear  moduli  in  the  1-3  and  2-3  directions,  respectively. 

If  nanal(2)=3  (sandwich)  then  the  material  is  a  sandwich  composite  and 
the  same  variables  entered  for  nanal(2)=0,2  are  entered  for  the  sandwich 
facesheets  only. 

•  CARD  11a:  els,  e2s,  gl2s,  pnul2s,gl3s,  gl2s  This  card  is  used  only 
when  nanal(2)=3.  The  above  properties  have  the  same  definition,  except 
they  apply  to  the  sandwich  core  material. 

•  CARD  lib:  If  nanal(2)=l,  (isotropic)  skip  this  card. 

If  nanal(2)=0,2  (composite),  read  in  np,  pt 

np  is  the  number  of  plies, 
pt  is  the  average  ply  thickness. 

If  nanal(2)=3  (sandwich  composite),  read  in  np,  pt,  pts 
np  is  the  number  of  layers  in  the  facesheets, 
pt  is  the  thickness  of  the  facesheet 
pts  is  the  thickness  of  the  core. 

•  CARD  11c: 

If  nanal(2)  =  l  (isotropic),  skip  this  card. 

If  nanal(2)=0,2,3  (composites),  read  in  the(i)  which  is  the  orientation  of 
the  plies  (in  degrees)  within  the  laminate.  For  the  sandwich  material,  it  is  the 
orientation  of  the  facesheet  plies. 

•  CARD  12:  If  iel=l  (flat  plate),  skip  this  card. 

Otherwise,  read  in  rad,  which  is  the  cylinder  shell  radius. 

•  CARD  13:  nfor  is  the  number  of  degrees  of  freedom  where  the  nodal  forces 
are  to  output. 

•  CARD  13a:  ifor(nfor)  is  an  array  of  nfor  global  degree  of  freedom  num¬ 
bers  where  the  nodal  forces  will  be  calculated  and  output.  For  example,  in 
displacement  control,  it  might  be  beneficial  to  know  the  resulting  forces  along 
the  edge  where  the  displacements  are  input. 

•  CARD  14a:  nstress,  ifail,  icrit 

nstress  is  the  number  of  elements  where  stress  calculations  are  to  take 
place.  If  failure  criteria  are  used,  this  number  will  be  automatically  (internally) 
set  to  the  number  of  elements. 
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ifail  determines  if  failure  analysis  is  used:  0  for  no  failure  analysis,  1  for 
failure  analysis. 

icrit  when  ifail=l,  this  variable  determines  the  type  of  criterion  used:  1 
for  Hashin,  2  for  Lee,  3  for  Maximum  stress. 

•  CARD  14b: 

If  ifail=l  (failure  analysis)  and  nanal(2)  =0,2,3  (composite  laminate)  read 
in  sigfn,  sigfs,  sigmt,  sigmc,  sigms,  sigdt,  sigds  corresponding  to  the 
material  strengths  of  the  chosen  composite.  The  fourth  letter  corresponds  to: 
f  for  fiber,  m  for  matrix,  and  d  for  delamination.  The  fifth  letter  represents:  n 
for  normal,  s  for  shear,  t  for  tension,  c  for  compression. 

•  CARD  14c: 

If  nanal(2)=3  (sandwich)  also  enter  the  core  material  strengths  for  failure 
evaluation.  The  required  strengths  are  sigcr,  sigl3c,  sig23c 

sigcr  is  the  transverse  direct  strength. 

sigl3c  is  the  transverse  core  shear  strength  in  the  1-2  direction. 
sig23c  is  the  core  shear  strength  in  the  1-3  direction. 
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1/4  of  24  x  24,  [0/90] 2s,  Tension,  Load  Control,  Use  Failure 
2,8,0,2,0,1,0,0,16 
0,11,80,0,0.001,1 
0.1,1.0,0.88 
12,12 

0.25,0.25,0.25,0.25,0.25,0.25,0.25,0.25, 

0.25,0.25,0.25,0.25,0.25,0.25,0.25,0.25, 

0.25,0.25,0.25,0.25,0.25,0.25,0.25,0.25, 

0.25,0.25,0.25,0.25,0.25,0.25,0.25,0.25, 

0.25,0.25,0.25,0.25,0.25,0.25,0.25,0.25, 

0.25,0.25,0.25,0.25,0.25,0.25,0.25,0.25 

105,106,107,108,117,118,119,120, 

129,130,131,132,141,142,143,144 

3,-17146. 

25 

1,26,39,64,77,102,115,140,153,178, 

191,216,229,254,267,292,305,330,343, 

368,381,406,419,444,457 

106 

25,1,0,0,1,0,1,0 

38,1,0,0,0,0,0,0 

63,1,0,0,1,0,1,0 

76,1,0,0,0,0,0,0 

101,1,0,0,1,0,1,0 

114,1,0,0,0,0,0,0 

139,1,0,0,1,0,1,0 

152,1,0,0,0,0,0,0 

177,1,0,0,1,0,1,0 

190,1,0,0,0,0,0,0 

215,1,0,0,1,0,1,0 

228,1,0,0,0,0,0,0 

253,1,0,0,1,0,1,0 

266,1,0,0,0,0,0,0 

291,1,0,0,1,0,1,0 

304,1,0,0,0,0,0,0 

329,1,0,0,1,0,1,0 

342,1,1,0,0,0,0,0 

367.1.1.1.1.1.1.1 
380,1,1,0,0,0,0,0 

405.1.1.1.1.1.1.1 
418,1,1,0,0,0,0,0 

443.1.1.1.1.1.1.1 
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456,1,1,0,0,0,0,0 
481,1,1,1,1,1,1,1 
458,0,1,0,0,0,0,0 
459,0,1,0,0,1,0,1 
460,0,1,0,0,0,0,0 
461,0,1,0,0,1,0,1 
462,0,1,0,0,0,0,0 
463,0,1,0,0,1,0,1 
464,0,1,0,0,0,0,0 
465,0,1,0,0,1,0,1 
466,0,1,0,0,0,0,0 
467,0,1,0,0,1,0,1 
468,0,1,0,0,0,0,0 
469,0,1,0,0,1,0,1 
470,0,1,0,0,0,0,0 
471,0,1,0,0,1,0,1 
472,0,1,0,0,0,0,0 
473,0,1,0,0,1,0,1 
1,0, 1,1, 1,1, 1,1 
26,0,1,0,0,0,0,0 
39,0,1,1,1,1,1,1 
64,0,1,0,0,0,0,0 
77,0,1,1,1,1,1,1 
102,0,1,0,0,0,0,0 
115,0,1,1,1,1,1,1 
140,0,1,0,0,0,0,0 
153,0,1,1,1,1,1,1 
178,0,1,0,0,0,0,0 
191,0,1,1,1,1,1,1 
216,0,1,0,0,0,0,0 
229,0,1,1,1,1,1,1 
254,0,1,0,0,0,0,0 
267,0,1,1,1,1,1,1 
292,0,1,0,0,0,0,0 
305,0,1,1,1,1,1,1 
330,0,1,0,0,0,0,0 
343,0,1,1,1,1,1,1 
368,0,1,0,0,0,0,0 
381,0,1,1,1,1,1,1 
406,0,1,0,0,0,0,0 
419,0,1,1,1,1,1,1 
444,0,1,0,0,0,0,0 


457,0,1,1,1,1,1,1 

339,1,1,0,0,0,0,0 

340,1,1,0,0,0,0,0 

341,1,1,0,0,0,0,0 

360,1,1,0,0,0,0,0 

361.1.1.1.1.1.1.1 
362,1,1,0,0,0,0,0 

363.1.1.1.1.1.1.1 
364,1,1,0,0,0,0,0 

365.1.1.1.1.1.1.1 
366,1,1,0,0,0,0,0 
377,1,1,0,0,0,0,0 
378,1,1,0,0,0,0,0 
379,1,1,0,0,0,0,0 
398,1,1,0,0,0,0,0 

399.1.1.1.1.1.1.1 
400,1,1,0,0,0,0,0 

401.1.1.1.1.1.1.1 
402,1,1,0,0,0,0,0 

403.1.1.1.1.1.1.1 
404,1,1,0,0,0,0,0 
415,1,1,0,0,0,0,0 
416,1,1,0,0,0,0,0 
417,1,1,0,0,0,0,0 
436,1,1,0,0,0,0,0 

437.1.1.1.1.1.1.1 
438,1,1,0,0,0,0,0 

439.1.1.1.1.1.1.1 
440,1,1,0,0,0,0,0 

441.1.1.1.1.1.1.1 
442,1,1,0,0,0,0,0 
453,1,1,0,0,0,0,0 
454,1,1,0,0,0,0,0 
455,1,1,0,0,0,0,0 
474,1,1,0,0,0,0,0 

475.1.1.1.1.1.1.1 
476,1,1,0,0,0,0,0 

477.1.1.1.1.1.1.1 
478,1,1,0,0,0,0,0 

479.1.1.1.1.1.1.1 
480,1,1,0,0,0,0,0 
0.,0.,0.,0,0.,0.,0.,0,0.,0., 
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0.,0.,0.,0.,0.,0.,0.,0.,0.,0., 

o.,o.,o.,o.,o.,o.,o.,o.,o.,o., 

o.,o.,o.,o.,o.,o.,o.,o.,o.,o., 

0.,0.,0.,0.,0.,0.,0.,0.,0.,0., 

0.,0.,0.,0.,0.,0.,0.,0.,0.,0., 

0.  ,0.,0.  ,0.,0.  ,0.  ,0.  ,0.  ,0.  ,0. , 
0.,0.,0.,0.,0.,0.,0.,0.,0.,0., 
0.,0.,0.,0.,0.,0.,0.,0.,0.,0., 
o.,o.,o.,o.,o.,o.,o.,o.,o.,o., 

0.,0.,0, 

o.,o.,o.,o.,o.,o.,o., 

o.,o.,o.,o.,o.,o.,o., 

0.,0.,0.,0.,0.,0.,0., 

0.,0.,0.,0.,0.,0.,0., 

o.,o.,o.,o.,o.,o.,o., 

0.,0.,0.,0.,0.,0.,0., 

0.,0.,0.,0.,0.,0.,0., 

0.,0.,0.,0.,0.,0.,0., 

0.,0.,0.,0.,0.,0.,0., 

0.,0.,0.,0.,0.,0.,0., 

0.,0.,0.,0.,0.,0.,0., 

0.,0.,0.,0.,0.,0.,0., 

0.,0.,0.,0.,0.,0., 

0.,0.,0.,0.,0.,0.,0.,0.,0.,0., 
0.,0.,0.,0.,0.,0.,0.,0.,0.,0., 
0.,0.,0.,0.,0.,0.,0.,0.,0.,0., 
0,0,0.,0,0,0.,0,0-,0,0., 
0.,0.,0.,0.,0.,0.,0.,0.,0.,0., 
0.,0.,0.,0.,0.,0.,0.,0.,0.,0., 
0.,0.,0.,0.,0.,0.,0.,0.,0.,0., 
0.,0.,0.,0.,0.,0.,0.,0.,0.,0., 
0.,0.,0.,0.,0.,0.,0.,0.,0.,0., 
0.,0.,0.,0.,0.,0.,0.,0.,0.,0., 
o.,o.,o.,o.,o.,o.,o.,o.,o.,o., 
o.,o.,o.,o.,o.,o.,o.,o.,o.,o., 
0.,0.,0.,0.,0.,0.,0.,0.,0.,0., 
0.,0.,0.,0.,0.,0.,0.,0.,0.,0.  0 
19.7e6,1.579e6,0.925e6,0.276,0.925e6,0.462e6 
8,0.00514 

0.,90.,0.,90.,90.,0.,90.,0. 

12.0 

25 
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1,116,142,257,283,398,424, 

539,565,680,706,821,847, 

962,988,1103,1129,1244,1270, 

1385,1411,1526,1552,1667,1693 

0,1,1 

216.1e3,31.7e3,6.31e3,6.31e3,6.31e3,7.20e3,12.39e3,17.39e3 
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14.  SUBJECT  TERMS 


composite,  cylinder,  shell,  cutout,  failure,  fatigue 


15.  NUMBER  OF  PAGES 
110 


16.  PRICE  CODE 


17.  SECURITY  CLASSIFICATION  18.  SECURITY  CLASSIFICATION  19.  SECURITY  CLASSIFICATION  20.  LIMITATION  OF  ABSTRAC 
OF  REPORT  OF  THIS  PAGE  OF  ABSTRACT 


UNCLASSIFIED 


UNCLASSIFIED 


UNCLASSIFIED 


Standard  Form  298  (Rev.  2-89  EG 

Prescribed  by  ANSI  Std.  239.18 

Designed  using  Perform  Pro,  WHS/DIOR,  Oct  94 


